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ABSTRACT 

An effective approach in finite element analysis of the stress 
field at the traction free boundary of a solid continuum was studied. 
Conventional displacement and Assumed Stress Hybrid finite elements 
were used in determination of stress concentrations around circular 
and elliptical holes. Specialized Hybrid elements were then developed 
to improve satisfaction of prescribed traction boundary conditions. 
Results of the stress analysis indicated that finite elements which 
exactly satisfy the free stress boundary conditions are most accurate 
and efficient in such problems. A general approach for Hybrid finite 
elements which incorporate traction free boundaries of arbitrary geom- 
etry was formulated. 
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CHAPTER 1 


Introduction 

The ''Finite Element Method" has proven to be an effective and con- 
venient approach to solution of problems of continuum mechanics . Orig- 
inally conceptualized in the 1950s as a tool in solving problems of 
elasticity and structural mechanics [1] , the finite element method has 
been increasingly used in many different fields of engineering/ applied 
mathematics, and physics. Rigorous mathematical formulation of the 
method based on different variational principles has afforded scientists 
greater flexibilities in the selection of the problem, formulation of 
the solution, and implementation of the method. Today finite elements 
are used in problems of heat transfer and fluit flow, as well as the 
traditional structural mechanics for which they were developed. 

The well known assumed displacement approach, based on the prin- 
ciple of minimum potential energy [2] , [3] is widely used in finite 

element formulations of structural mechanics. Although a successful 
method in many applications, inherent drawbacks of the formulation do 
exist. In the assumed displacement formulation the requirement of dis- 
placement compatibility is identically satisfied throughout the element, 
by selection of an appropriate set of displacement interpolations. This 
selection process is not always guaranteed and in fact has posed dif- 
ficulties in construction of conventional plate bending elements. In 
the displacement model, equilibrium of applied forces and the resultant 
stresses is in general not satisfied except in the limiting case where 
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more and more degrees of freedom are used in modeling the structure 
until monotonic convergence is achieved with compatible elements. Since 
in the assumed displacement formulation stresses are obtained by dif- 
ferentiation of the displacements, any error of the displacement solu- 
tion is amplified in this process and will hence introduce distortion 
or noise in the stress or strain solution. It is for this reason that 
when stress concentrations exist in the vicinity of structural discon- 
tinuities or supports, a finer mesh of elements is always used. The 
compatible assumed displacement approach will in principle underestimate 
the displacement solution by overestimating the stiffness of the struc- 
ture . 

On the other hand, finite elements can be formulated based on the 
principle of minimum complementary energy [3] * In this approach the 
stresses satisfy the equilibrium conditions and become the unknown 
variables directly solved for. In principle the minimum complementary 
energy formulation will overestimate the displacement solution by under- 
estimating the structural stiffness. Since the condition of displace- 
ment continuity along interelement boundaries is not satisfied, the 
formulation will generate incompatible finite elements. 

As introduced by Pian [4] [5] , the assumed stress hybrid formula- 
tion relaxes displacement interpolation requirements critical to the 
assumed displacement model, while assuring the interelement compat- 
ibility condition not achieved by the equilibrium model. Rigorously 
formulated by Tong and Pian [6] , the hybrid model is based on a modified 
complementary energy principle in which stresses are the variables 
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directly solved for. 


1.1 Objective 

For an accurate estimation of the stress field at the streee free 
boundary of a solid continuum, it is desirable to choose finite elements 
that satisfy the prescribed stress boundary conditions. For such prob- 
lems, hybrid stress formulation becomes a more logical alternative be- 
cause requirements of stress equilibrium and free boundary traction can 
be simultaneously satisfied. 

To assess the effectiveness of assumed stress hybrid elements under 
prescribed free traction conditions, stress concentrations around cir- 
cular and elliptical holes subjected to unidirectional tension were 
examined. For the purpose of this investigation, three distinct hybrid 
elements are studied. The first study examines the effectiveness of 
regular isoparametric hybrid elements that do not meet the condition of 
free boundary traction. The performance of these elements are then 
judged against conventional displacement based isoparametric elements 
with regard to solution accuracy and efficiency • Next a second class 
of hybrid elements that approximately satisfy the prescribed free- 
traction boundary condition, via the boundary point matching technique, 
is considered. Finally, the advantage of using hybrid elements that 
identically satisfy the traction-free boundary condition is examined 
in the special case of a circular boundary. In each case finite ele 
ment stress solutions are obtained for a comparative study against the 
analytical solution. A general formulation for assumed stress hybrid 
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finite elements that incorporate a traction-free boundary of arbitrary 
geometry is presented. 
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CHAPTER 2 


The Assumed Stress Hybrid Element 
Formulation of the assumed stress hybrid element is in general 
based on either the stationary value of the modified complementary 
energy 7 T MC , or of the Hellinger-Reissner principle [3]. It can be 
shown that for a given set of assumed stresses these formulations pro- 
duce identical results [5]. Thus depending on the problem at hand, 
ease of implementation or economical considerations may prove one 
formulation advantageous to the other . 

2.1 The Modified Complementary Energy Principle 

In the modified complementary energy principle, stationary value 
of the functional 77 ^ is established [3]; 

/*T C 

7 L r . z u f <r T srdv. [r T uds J f T u ds l 12.11 
” \ z \ ~ ~ ~ Jw - Js '" - 1 

where: 

* Number of elements. 

Vn s Spatial domain of element n. 

^^5 Boundary of Vrj . 

- Portion of 314 where tractions are prescribed. 

(/ = Self equilibrating stresses. 

^ x Compliance matrix. 

7* = Traction matrix. 
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£5 £ Displacement matrix. 

T i Prescribed tractions. 

Although in the general formulation body forces are also accounted for, 
their presence is intentionally neglected in equation (2.1). The 
existence of body forces will be considered in the Reissner formulation 
in a much simpler manner identical to that of the assumed displacement 
model . 

In the Finite Element formulation the stresses (T are assumed over 
the entire domain, and the displacements are interpolated from the 
nodal values along the boundaries of the element only. 


f * 

f * 


where : 


p 




/ 


Matrix of assumed polynomials. 
Undetermined coefficients. 
Interpolation functions. 

Nodal displacements . 


In Equation (2.1) the product of the traction vector and the displace- 
ment vector^ is integrated along the boundary of the element. On any 


6 


boundary the traction vector can be determined directly from the assumed 
stresses (f \ 


/ - y T a * v r P fi 


(2.4) 


where : 

y - Matrix of direction cosines. 

A' 

Substitution of equations (2.2), (2.3) and (2.4) into equation (2.1) 

will discretize the functional « 


n | z L 


(2.5) 


y 

v> 


/ 


-A 7 / - f j 

/ A J S<r J- - 1 

v /7 


Redefining the integrals of equation (2.5); 


H r /" P r s p dv 

Jtfn ~ ~ ~ 

6 = f p 7 y L 

hv* 



(2.6) 


(2.7) 


( 2 . 8 ) 
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2.2 The Hellinqer-Reissner Principle 


Assuming satisfaction of displacement boundary conditions, the 
Hellinger-Reissner principle requires the stationary value of the 
functional [3]; 


v* - % |-jf [/if - f 


(2.9) 


where j is now the matrix of prescribed body forces. Here the inte- 
gration of the stresses is performed in the spatial domain and matrices 
& and 5 are identical to those defined for equation (2.1). The 
essential difference between the and / 7j^ approach is in the inter- 
polation of displacements. Where in the '77 ^ approach displacements 
are interpolated only along the boundaries (2.3), in the '7^ formula- 
tion they are interpolated over the entire domain; 


a * 



( 2 . 10 ) 


and matrix & is a differential operator which by operating on the dis- 
placements gives the strain matrix; 


D u 



( 2 . 11 ) 


Substitution of equations (2.2) and (2.11) into equation (2.9) will 
discretize the functional 7Tg ; 
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( 2 . 12 ) 



Redefining the last three integrals of equation (2.12) 

<3 r / dV (2.13) 

J\Zn 

Q 7 f / A dV - f 7 7 A dS 

J V»~ ~ J S<r~ ' 

'v* > 1L * ^ T gy -Q r y j < 2 - i4 > 

Hence equivalence of the two functionals (2.8) and (2.14) is apparent. 
By virtue of their mathematical equivalence, body forces present in 
the 'Tffic formulation can be accounted for in matrix Q, exactly analogous 
to that of the ^ formulation [5] . 

2.3 Solution of the Discretized Equations 

In both the modified complementary Energy and the Hellinger- 
Reissner formulation, the functional was seen to reduce to the form; 



Z ,i r », j ; Mi ~ Q r i 1 

* ) /- / i ~ L ) 


(2.15) 


9 


The functional has to assume a stationary value with respect to the 
independent variable ; 


'dTTu 

% 


- o 



( 2 . 16 ) 


Substitution of equation (2.16) into (2.15) gives; 

Variation of ^ with respect to ^ gives; 

( § a ~ ) ? * ® (2-i 

Comparison of equation (2.17) with the familiar assumed displacement 
results; 

h - r 


establishes the equivalent stiffness matrix obtained by the Assumed 
Stress Hybrid method 



( 2 . 18 ) 
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Matrix (5 of the above equation, depending on the formulation used, is 
defined by either equation (2,7) or equation (2.13). 

2.4 Constraint Conditions on the Assumed Stresses 

The assumed stresses as defined by equation (2.2) are not entirely 
arbitrary and have to satisfy some minimum requirements. 

1) The assumed stresses have to satisfy the homogeneous equa- 
tions of equilibrium, that is the equations of equilibrium 
for no applied forces. In most problems of continuum mech- 
anics the satisfaction of the equilibrium conditions can be 
achieved by extraction of the stresses from an assumed stress 
function. 

2) Another condition imposed on the assumed stresses is the 
minimum number of betas required to avoid introduction of 
spurious energy modes. If N is the number of degrees of 
freedom present in the element, and M the number of 
rigid body motions the element can undergo, then at least 
N-M betas have to be used in the formulation of the stiffness 
matrix. 

2.5 Computation of the Stresses 

Stresses can be calculated directly from equation (2.2); 


<r, p i 
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where 



is given by equation 


(2.16) 


Jk - f 

O' - P H 6 j (2.i9) 

The matrix product H 6 will be computed during the element stiffness 
computation (equation 2.18), and can be stored for the stress calcula- 
tion, after the displacement solution is obtained. 

2 . 6 Implementation Considerations 

It is evident from equation (2.6) that the size of the matrix H is 
directly related to the number of betas used in the assumed stress 
matrix; for example, the size of matrix H will be n x n if n betas are 
utilized. Since matrix N has to be inverted (equation 2.18) for the 
determination of the element stiffness matrix, it is desirable to keep 
the number of betas to a minimum for computational efficiency. However, 
when polynomials are used in the expression of the assumed stresses, a 
complete polynomial is in general more desirable than an incomplete 
expression randomly truncated for efficiency sake. This is due to the 
fact that a complete polynomial will render an invariant element stiff- 
ness matrix, while an incomplete polynomial will obviously affect the 
element stiffness matrix based on its orientation relative to the 
Global system of coordinates. This problem is of course alleviated 
if a local system of coordinates is used for the stiffness computation, 
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although a transformation of the stiffness matrix from the local to the 
global system of coordinates would then be necessary. In plane elas- 
ticity problems it is possible to impose further restrictions on the 
assumed stresses so that they also satisfy the compatibility conditions 
throughout the element [7] . The constraints reduce the number of betas 
used in a set of complete polynomials and maintain element invariance 
at the same time. This approach is also deemed desirable since the 
element will then satisfy both the condition of stress equilibrium and 
compatibility, but is in general not easily applicable especially in 
three dimensional solid elements. 

The difference between and ^ based hybrid element amounts 

to the difference in the computation of the matrix (5 (equation 2.7 and 
2.13). Since equations (2.7) and (2.13) are in general integrated 
numerically, it would appear that the model would be more efficient 
since the integration is only performed around the boundaries of the 
element. However, the comparison also requires an investigation of the 
order of terras being integrated. In the 'TT# formulation the matrix 3 
(equation 2.13) is obtained by differentiating the displacements 
(equation 2.11) and therefore the product P 3 is lower in order than 
the product for the same matrix P. The most efficient approach 

will thus depend on the number of degrees of freedom of the element and 
the highest order of powers appearing in the assumed stresses. Another 
consideration may arise from the construction of the matrix of direc- 
tion cosines V . In some general three dimensional elements with 
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boundaries of arbitrary geometry it may be very difficult to establish 
V while this problem is nonexistent in the 17% formulation. 

Finally, it should be noted that the only rigorous requirement on 
the displacement interpolation matrix A of equation (2.10) is that it 
reduce to the proper interpolation of the displacements along the 
boundaries of the element, and is otherwise arbitrary. This could 
somewhat relax the determination of the interpolation functions for an 
element with uneven distribution of nodes. 
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CHAPTER 3 


Example Problem 

To determine the performance of Assumed Stress Hybrid finite ele- 
ments a plate of finite width with a circular or elliptical hole is 
analyzed. The plate is subjected to uniaxial tension and hence a con- 
dition of stress concentration is developed at the hole, across from 
the center of the hole directly perpendicular to the direction of the 
applied load. Analytically obtained values of the stress concentra- 
tion factor are known from the plane elasticity solutions with both 
circular and elliptical holes [8] . 

Figure 1 shows the dimensions of the plate with a circular hole 
subjected to tensile stress^? . A stress concentration factor of 4.32 
exists at points A and B on the circular hole. Figure 2 shows the 
dimensions of the plate with the elliptical hole subjected to the same 

tensile stress p . For this configuration a stress concentration factor 

r o 

of 9.50 is present at points A and B shown on the ellipse. The length 
of the plates compared to the dimensions of the holes were judged long 
enough to simulate a uniform tensile pressure acting at infinity. 

The two dimensional continuum nature of the problem selected for 
analysis is important since the results will have a direct bearing on 
similar analyses done with three dimensional continuum elements. The 
solid three dimensional continuum element can in fact be regarded as 
a natural extension of its two dimensional continuum counterpart. 
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Circle; x 2 +y 2 =1.0 

Figure-1 Plate With A Circular Hole 



Ellipse; x 2 +(^j) 2 =1.0 


Figure-2 Plate With An Elliptical Hole 




Due to the symmetric nature of the problem only a quarter of the 
plates were used for analysis, namely the shaded regions of Figures 1 
and 2. On the boundaries where the plates were artificially cut, zero 
displacements were prescribed in directions perpendicular to each 
respective boundary. 

For the finite Element Analysis of the plates, the coarse mesh of 
Figure 3, labelled Mesh-1, was first used. In Mesh-1 the hole is 
divided into two sections by a line from the center of either the circle 
or the ellipse intersecting at 45 degrees. Stresses were computed at 
locations A shown in Figures 3a and 3b, and were compared against the 
expected results. It should be noted that Mesh-1 incorporates a total 
of three elements and in general is not expected to render accurate 
results. 

A finer mesh can be obtained by subdividing Mesh-1. This is 
achieved by connecting the mid-boundaries of all the elements in 
Mesh-1, which in effect replaces each element with four subelements. 

This mesh is labelled Mesh-2 and is shown in Figure 4. Mesh-2 incor- 
porates a total of twelve elements and should in principle yield a 
better solution compared to Mesh-1. 

The Finite Element Analysis Basic Library (FEABL) of the Aero- 
elastic and Structures Research Laboratory of the Massachusetts Insti- 
tute of Technology [9] was used for assemblage of the elements and 
solution of the equations. 
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Figure-3a Mesh-1 ; Quarter Plate With 
Circular Hole 



Figure-3b Mesh-1 ; Quarter Plate With 
Elliptical Hole - 
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Figure-4b Mesh-2 ; Quarter Plate With 
Elliptical Hole 
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CHAPTER 4 


A 4-8 Variable Node Hybrid 
Isoparametric Quadrilateral Element 


4. 1 Objective 

The displacement based eight node isoparametric element has proven 
very effective in two dimensional continuum problems such as plane 
stress or plane strain. The availability of the eight node plane 
element in most of today's advanced Finite Element Library codes con- 
firms its extensive acceptability. To this end it is therefore logical 
that for analysis of the problems described in Chapter 3, with Finite 
Element codes utilized by the industry, the eight node displacement 
based element would be the effective option. 

Development of an eight node hybrid isoparametric element will 
provide a rational alternative for comparison against the eight node 
displacement based element. The Hybrid Isoparametric Quadrilateral 
element developed (called HISQUE) is formulated such that any number 
of the mid-nodes can be removed. The number of stress terms (i.e., 
the number of betas) programmed in the element is also user defined 
as an added flexibility for an optimum selection based on the number 
of nodes employed. The element HISQUE can employ up to a complete 
cubic set of assumed stresses in the XY plane. 
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4. 2 Formulation of Element HISQUE 


4,2.1 The Assumed Stresses 

Element HISQUE is formulated based on the Hellinger-Reissner 
principle: 


M = [ P r S P dv 
0 - f P r 3 dv 


The matrix of assumed stresses is defined: 

(Tx* 

(fa > 

(7xy 

and can be determined from the Airy stress function 



( 2 . 6 ) 


(2.13) 


(4.1) 


<r» - 

'di 1 

cr n . 

<ui , - 

a*?y 


(4.2a) 


(4.2b) 


(4.2c) 


As defined by equations (4.2), the stresses obtained will satisfy the 
homogeneous equations of equilibrium in Cartesian coordinates: 

9jT*y i Wyy = o 


7* 7Y 
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o 


'd < Tw g£>y 
'd* f dy 

If an infinite series of polynomials in ascending powers of x and 
y is used for the stress function $ , where the individual terms of the 
series are weighted by arbitrary constants (Xi 

cf> =■ cx o + a, x -f «z y * « 3 X ; f a .4 xy + cc s ... (4.3) 

then substitution of equation (4.3) into equations (4.2) will determine 
the stress vector (4.1). The coefficients ^ are then renamed y^,' such 
that t corresponds to the lowest coefficient OC^ appearing in the 
series . Since the element HISQUE can have a maximum number of eight 
nodes with two degrees of freedom at each node, and also since the 
element can undergo three rigid body motions, it follows that a mini- 
mum of 13 terms are required in the assumed stresses (Chapter 2) . The 
stresses obtained from the Airy stress function (4,3), 

(Txt s fit + y ji 4 + + yyl 6 + (4.4a) 

/ 2xy y^U 1 - /y^/j + Z 2 y^is V xy 2 yiji 

t J + 
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(4.4b) 


fay c /h + x^6 4 y jU + x 2 ^ „ y^ /0 
+ 2*y fit. z + X 2 /fa 4 + 3x y 2 /l,s + L 3 
4 X 2 y ^ /S 4 '•■ 


fay * ^>3 _ y fit . //1 7 _ ZXY^JO - 

- x * ft i 2 - 3x 2 y j^> ts - y/^ " (4 * 4c) 

" 5 ^ /^ /l3 y " ' 

or in matrix form, 


fax 


r 


i 

- 

i 

< fay 

> - 

P 

< 

: f 

fay 

. 4 


- * 


/^/fl 
> . 


From equations (4.4) it is shown that the highest power appearing in 
the stresses depends on the number of betas retained. As seen in 
equations (4.4) , using eighteen betas will result in a complete cubic 
expression of the stresses. 

4.2.2 The Strain Matrix 

The strain matrix 3 of equation (2.13) is obtained by the dif- 
ferentiation of the displacements 
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a - 


(4.5) 



£ s 



(4.6) 


The matrix ,4 of equation (4.5) is the matrix of the interpolation func- 
tions. In the isoparametric transformation of the coordinates the fol- 
lowing relations hold: 


/ * Z h 

/C 

y* ? - ^ 

where : 

h^ * Interpolation functions. 

(x.,y.) a Nodal coordinates. 

1 x 

or in matrix notation: 

Xi 

y» 



(4.7a) 


(4.7b) 
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Since in the isoparametric transformation the displacements are inter- 


polated identically: 

U x , X- h; a * 

U y s 71 h * tfi 

i 

where u^ and are nodal displacements of node i in x and y directions 
respectively. It follows directly that matrix //of equation (4.7c) is 
identical to matrix 4 of equation (4.5). 

For the eight node isoparametric element shown in Figure 5, the 
interpolation functions h^ are [10] : 

h, * l (Ui)(UD-l h 5 h s - j. 
h z -- l(i- l)('+'i)-± 6s6s - 

h 3 * l('-l)(i-y)-l S 6 h(. _ L$>?h 7 

h<,L(Ut)0-l) .jiaht M . B) 

hi: Sg.(l-l l )(Ul) 
hi : k(>- l)(l-q l ) 

h? & 0-12)0-1) 

- S*(U DO- j z ) 

^ I / J f nodi jt pmt/ir 

( 0 // node -c acse./?^ 

By employing the on/off switches & , any combination of the nodes 
numbered 5, 6, 7 and 8 can be removed (Figure 5) while an appropriate 
displacement interpolation is preserved. 
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The interpolation functions are functions of £ and ^ which in 
turn means that the displacements as expressed by equation (4.5) are 
also functions of £ and-^. Derivatives of the displacements with 
respect to the x, y coordinates can be determined as follows: 


dit; _ 2a. 7^ . 9a> 2L 

Vt " 2* 71 V 

W „ To, 2a. 21 

° r ' 


' 9 ‘ 

91 

> a: = 

' j 

< 

9 ' 
> 

2 

n t 


*** 


! ?n 

j 


where J is the two by two Jacobian matrix, the entries of which can 
be easily obtained by differentiation of equations (4.7a) and 4.7b). 
From equation (4.9) it follows: 


* 7 ' 

< 


‘ J ' 

_ 1 

< 

^ » 

_ j 


*** 

„ J 


2. 


where : 


J 

_ 1 

1 

Jzz 

-J,Z ‘ 


1 J 

i 1 ~ 

_ J/2 

J// 


The B matrix can then be written as: 
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dht 

0 

d/n 

0 

* • « 

0 




9% 




3 = 


dht 

0 

%2 

• # • 

dhs 



djj 




¥ 


2hi 

'dht 

9hi 

9 hi 


9ha 


[ W 

21 

% 

21 


dx 

where 

[from 

t equation (4.10] 

I] : 





n 

l * 

J_ ( Jit iJlL 

171 ' dl 

- J/z 

dhi \ 
drj 1 



dA, 

— (- 

Jjl 'ihi . 

+ J/t 

Vhi) 



dy 

131 l 

21 


97f / 


4.2.3 Numerical Integration and Implementation 

The matrices /V and 6 of equations (2.6) and (2.13) are 
by the Gauss-Quadrature numerical integration scheme [11] : 


7 , 



f(X,Y) dv 



d A 


A 


t r Element thickness; assumed constant. 
dA = Differential area of the element. 

Equation (4.11) becomes: 

+ ! + / 



(4.11) 


determined 


(4.11) 


(4.12) 
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where : 


^/nUlsn *r$jr?dr:± 

D 


(t • / 1j ) S 6 ju£C lnh ; ^ 

W‘ ) Ml Z &,3-J££ wtijhhf!/ f"X chr$ 

/ o c / 


Using n points in the Gauss integration scheme of equation (4.12) 
(i.e., M = N = n) , a polynomial of order 2n - 1 can be integrated 
exactly. The highest power of x or y that appears in equation (2.6) 
will depend on the number of betas used in the stress terms. If the 
complete cubic expression of equations (4.4) with eighteen betas is 
used, the product P r $ P will contain polynomial terms of sixth order 
and hence a 4x4 Gauss quadrature scheme would seem adequate. For an 
accurate integration of the product P*3 , the number of nodes need to 
be considered. Equations (4.8) show that the interpolation functions 
h^ are in general second order in ^ and ^ . From equation (4.11) it 

is seen that typical terms appearing in & are of the form: 



where S replaces either^ or/? . 

When all nodes are present, the matrix will contain fourth powers 
of £ and/or//. With a complete cubic distribution of the stresses it 
is again seen that a four by four Gauss integration will be sufficient 
to accurately integrate ^7. Naturally the required number of integra- 
tion points will change if a lower order of stresses is used or the 
mid-nodes of the element are completely eliminated. 

Matrices M and 5 are hence determined as follows: 
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(4.13) 


H . ZZ [(f^JI/IL. W: VJj 
o , ZZ [(f 7 ^ M W/ 


In equation (4.13), S is the compliance matrix of an isotropic homo- 
geneous material in plane stress: 

' / - y o' 

s * J_ 

£ -v 1 o 

o o z(t+y) 

In equation (4.14), the determinant of the Jacobian matrix J resulting 
from the coordinate transformation is eliminated, because as seen in 
equation (4.11) the matrix has the reciprocal of this determinant 
as a scalar multiplier. 

The rationalization of the required number of integration points 
was arrived at experimentally. It was observed that, for example, 
using a5x5or7x7 Gauss integration scheme, makes insignificant 
changes to the stiffness matrix of the eight node element with com- 
plete cubic stress distribution. This deduction is substantiated in 
a report by Spilker [7] for an eight node hybrid isoparametric element. 
It should be noted that for strictly two dimensional analysis, imposing 
the condition of stress compatibility: 



will not only improve the efficiency of the element by reducing the num- 
ber of betas, but will also produce a stiffness matrix that models the 
physical problem more accurately [7]. However, in three dimensional 
analysis the conditions of stress compatibility do not appear as simple 
constraint conditions on the stresses and in general are very difficult 
to meet. For this reason, performance of the eight node element with 
the assumed stresses of equations (4.4) would be a better judge of the 
solid isoparametric element. 

The stiffness matrix of the element is then determined from equa- 
tion (2. 18) : 

K , 6 t h 

Since the inversion of matrix H can add considerable computation 
time, the most efficient approach should be considered. Many efficient 
routines are developed that factorize the matrix into triangular matrices 
and then invert it by backward and forward substitution. For element 
HISQUE, the subroutine DSINV, of the Scientific Subroutine Package (SSP) 
programmed by IBM, was used to invert H. This subroutine operates in 
double precision accuracy and requires storage for only the upper half 
of the symmetric matrix//, for inverting it by the Cholesky Factonza- 

tion method. 

The stiffness matrix was observed to have three zero eigenvalues 
corresponding to the three rigid body motions and hence free of spurious 

energy mode . 
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4. 3 Results and Discussion 
4.3.1 Plate with Circular Hole 

Stress solutions were obtained with a displacement based eight 
node isoparametric element (ELEM8) [12 ] at the point of stress concentra- 
tion with both Mesh-1 and Mesh-2. Satisfactory stresses (+5%) were ob- 
tained in both cases, compared to the expected solution of 4.32 at the 
point of stress concentration. The results are shown in Figure 6, where 
modeled with eight node elements Mesh-1 incorporates 36 degrees of 
freedom, and Mesh-2 106. The precedure was repeated with element HISQUE 
and the results are also shown in Figure 6. The error in the stresses 
is summarized in the table below: 



An alternative to modeling the problem entirely with displacement 
or hybrid elements is to use the hybrid element only where the stresses 
are critical. This treatment of hybrid elements as specialized elements 
is important from both solution efficiency and dependability aspects. 
Since the hybrid approach does not differentiate the displacements for 
the stress solution, it can be argued that where stress determination 
is the paramount objective, its utilization will become advantageous. 

A solution was obtained with Mesh-2, with element HISQUE used in 
the shaded section of Figure 4a, and element ELEM8 in the remaining part 
of the quarter plate. The percent error in the maximum stress for this 
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Figure-6 Maximum Stress (Lt In Plate With 


A Circular Hole 


case was -.9%. The difference in the error between this and the two 
previous cases, where complete hybrid and complete displacement elements 
were used, appears too academic. 

Indeed, except for the appearance of monotonic convergence of the 
hybrid stress solution versus the oscillating displacement based stress 
solution, Figure 6, use of the hybrid element HISQUE does not seem 
justified in the case of the circular hole. Nevertheless, with the finer 
mesh, better than 2% improvement is observed in the stress solution when 
element HISQUE is used. Closer examination of the stress equilibrium 
at the left boundary of the quarter plate. Figure 7, reveals that the 
resultant stress on this section is closer to the applied force when a 
complete hybrid, or partially hybrid - partially displacement based model 
is used. With Mesh-2 the following values were obtained for the 
resultant stress N: 


N s 1 (T xx dy 
-'s 

(4.16) 

S 3 Ooundartj of pl(XT£. 

A / D = 0.961 ~p 

(ALL EJL&M&) 

- /• 00 Q f o 

( A L.L rL/SCL’JL ) 


(ml asp ,v:c£jl; 


From Figure 7, it is also seen that, at the boundary, the stress 

solutions of the complete hybrid and the mixed model are indistinguish- 
able. 
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Figure-7 Stress Profile At The Left Edge With 
Eisht-Node Isonarametric Elements 


4.3.2 Plate with Elliptical Hole 

Stress solutions were obtained at the point of maximum stress of 

the plate with elliptical hole shown in Figure 2. Modeled according to 

Figures 2b and 4b, labeled Mesh— 1 and Mesh— 2, respectively, elements 

HISQUE and ELEM8 were both used in analysis. 

Analysis with the displacement based element ELEM8 did not produce 

acceptable results (+5%) with either mesh. As noted earlier, Mesh-1 

is really too coarse to realistically model the problem within finite 

element approximation, and its inclusion is merely intended to provide a 

convergence trend. With Mesh— 1, ELEM8 predicted a stress concentration 

factor of 7 . 29^p , an underestimation with more that 23% error compared 

to the expected value of 9.50 p . with Mesh-2, a stress concentration 

' 0 

of 8.57 was found which carries an error of -9.8%. 
o 

Analysis with element HISQUE converged with Mesh-2 and the results 
are summarized in the following table: 


Element 

error 

Type 

Mesh-1 

Mesh-2j 

ELEM8 

-23-3 

-9.78 | 

HISQUE 

-41.9 

-4.52 


The stresses obtained are shown in Figure 8. 

As done in the case of the circular hole, this plate was also 
analyzed with a combination of the two elements. Here, element HISQUE 
was used in the shaded region of Figure 4b, and the rest of the plate 
was modeled with displacement based elements ELEM8. A stress concentra- 
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factor of 9*04 was obtained in this analysis, which only carries an 
error of 4.84%. 

It is seen that employment of the hybrid element provides consider- 
able improvement over the strict displacement based element in the anal- 
ysis of the elliptical hole problem. Considering that the next refine- 
ment of the model would require 48 elements, the introduction of only 
five hybrid elements around the hole and on the left boundary becomes 
justifiable. So long as the cost of five hybrid elements equals that of 
41 displacement based elements, the reduced burden of modeling the plate 
with a fine mesh alone would rationalize the proposition. A more accu- 
rate comparison should also include the resulting increase in the global 
s t^ff ness and band-width, and hence the additional storage and solution 
cost. Although the ratio of required CPU (Central Processing Unit) time 
for analysis with elements HISQUE and ELEM8 is well below 8:1, there is 
much room in both for further improvements and a numerical comparison 
would not be a creditable representative. 

Figure 9 shows a plot of the stress distribution in Mesh- 2, along 
the left boundary, for all three cases. Again the indistinguishable 
results between the complete hybrid model and the mixed hybrid displace- 
ment model is emphasized. The reaction equilibrium can be determined 
from the resultant stress N of equation (4.16): 

/. O! -p ( ALL E~£MQ') 
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The fact that on the left edge, equilibrium is almost exactly satis 
fied and the stress distribution is nearly continuous across the two 
elements may erroneously suggest that the displacement solution of the 
stress field has attained convergence. 

A close examination of the normal stresses around the hole will 
P rov i<i6 an explanation of the more accurate hybrid analysis. The stress 
normal to the elliptical hole can be determined by transforming stress 
components along directions normal and tangent to the ellipse. The 
angle of the normal Of , is given by: 


<X 


/do 0 + 



tan 0 

. Os ) 2 


where : 

0 2 Angle of point (x, y) on the ellipse, in polar coordinates. 

Since the hole is not subjected to any external loads, a zero normal 
stress condition is in effect prescribed at the hole. Figure 10a shows 
a plot of the normal stresses obtained with the hybrid and the displace- 
analysis. The results obtained with the mixed model are again in- 
distinguishable from the complete hybrid analysis and hence are not shown 
^ evident from Figure 10a that the hybrid element approximates the 
stress boundary condition more accurately than the displacement element. 
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Figure-10a Distribution Of The Normal 

Stress At The Elliptical Hole 
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In terms of the absolute value of the stress deviation from zero, the 
displacement based results are at best 21% worse than the hybrid result, 
and over 700% at 0 = 22.5°. Figure 10b shows the normal stresses to the 
elliptical boundary normalized by the hoop stress. 

4 . 4 Concluding Remarks 

The ability of the hybrid element HISQUE to better satisfy stress 
boundary conditions around a hole, compared to the displacement based 
element ELEM8, is seen to produce better stress results in the cases of 
circular and elliptical holes. 

The most dependable, accurate and economical model in the finite 
element analysis of these problems is judged to be one that employs 
hybrid elements where accurate stress predictions are critical. The 
mixed model employs displacement based elements for the rest of the 
structure . 
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CHAPTER 5 


A- 7- 10 Variable Node Hybrid 
Isoparametric Quadrilateral Element 

5.1 Introduction and Objective 

Comparison of results obtained with eight node isoparametric ele- 
ments indicates that satisfaction of the free traction boundary condi- 
tion is essential in successful stress concentration analysis around 
holes. For this purpose, the formulation of a specialized hybrid iso- 
parametric element, HEL710, based on the collocation method suggested 
by Atluri [13], is presented. This approach led to introduction of two 
additional nodes on one boundary of the element HISQUE. For preserva- 
tion of compatibility in analysis with other elements the mid- nodes of 
the other three edges of the element are removable. The optimum stress 
distribution for this element was found to be a complete quartic. 

The objective is improvement of solution accuracy and efficiency 
compared to that obtained with a mixed model employing hybrid and dis- 
placement based eight node isoparametric elements. Element HEL710 is 
shown in Figure 11 in the Cartesian global coordinates. 

5.2 The Boundary Point Matching Technique 

As developed by Atluri and Rhee [13], the boundary point matching 
technique works by introducing the Lagrange multiplier u in the varia- 
tional formulation to enforce traction equilibrium. The Lagrange multi- 
plier u is then the compatible interelement boundary displacement. 
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n 



(5.1) 


frt? dV * f TudS+ f (T.T) % dS 

x -V JdY Js* 


~ The hybrid stress functional. 
ti 2 Compatible interelement boundary displacements. 

Portions of boundary where tractions are prescribed. 

T 2 Prescribed tractions. 

Comparison with equation (2.1) shows that the last integral is added to 
satisfy the traction equilibrium. If the following assumptions are 
made : 


O’ s P & (2.2) 

- - i: f < 2 . 3 ) 

U s P (K (5.2) 

Then one can write the traction vector in terms of the stresses: 


7,^'S y T f/i 


In equation (5.2) X are undetermined coefficients and/- are arbitrary 
order polynomials. 
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Substitution of equations (2.2), (2.3), (2.4) and (5.2) into 


equation (5.1) gives: 


r » - - - - - 


(5.3) 


The 


functional 77 ^ should attain a stationary value with respect to the 


two independent variables X and/4 . 


H 




C T oc ~ o 


(5.4a) 


*4 - 9 * 


(5.4b) 


where : 


Q r / F T ? dS 


and £ is determined by substituting coordinates of points on the boundary 
in the/ 5 matrix. From equation (5.4a): 


_ / . i ^ -r 

H C oc 


A r H 6 3 , h 


(5.5) 


Therefore, from equation (5.4b) : 
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9 r £H~'$ f t £ T * 


~ = - (££'*£ r ) ~ f + (cH m/ c r ) ~*q 


(5.6) 


Substituting (5.6) into (5.5): 


/i f‘* f - tt '£■ ?f r f r ; ' / 

By substitution of (5.6) and (5.7) into (5.3) and comparison with the 
assumed displacement functional 


% - ± 2 r tl 


£ J ~ / 


^i^roent stiffness matrix is obtained: 


i - SY's - 


(5.8) 


However, if an element with n nodes is considered at m nodes of which 

(n > m) zero traction is prescribed, then the general stiffness matrix 
°f this element 


48 



(2.18) 


^ 

can be rewritten as: 



7 1 

6a 


6 a 


*** 

H 


* 

1 

k ^0 

1 


6g 


Where & 3 is that portion of 6 corresponding to the m nodes at which 
tractions are prescribed to be zero. From equation (5.9) the stiffness 
matrix is seen to partition accordingly: 




Here>£,.. is a matrix of order (m x m) . 

1+0 O 

nodes B then: 


] 

I 

/<■ 83 : 


J 

Since no loads are applied at 
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Static condensation of the ^ degrees of freedom gives: 

?£ 


kr 3 
~ 7 A 
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where : 



- t! £a - !? a H ^3 (&6 H 6&)~ Gt H~ Ga 


(5.10) 


The equivalence of equations (5.8) and (5.10) is seen when 

£ ' i 

In other words the method of boundary point matching is equivalent to 
introducing additional nodes on a traction free edge and then statically 
condensing these nodes in obtaining a reduced stiffness matrix. 

Although very flexible, the boundary point matching technique has 
inherent drawbacks. The major cost- related criticism directed at the 
assumed strees hybrid formulation is the necessity of a matrix inversion 
for stiffness matrix evaluation. As is seen from equation (5.8) the 
boundary point matching requires a second matrix inversion, namely of the 
matrix product (CHC ). Furthermore, the size of this matrix is 
directly related to the number of collocation points used to satisfy the 
boundary traction condition; hence, the more accurate a solution desired, 
the higher a price for the element. 

To avoid costly inversions, element HEL710 employs only two addi- 
tional nodes on one boundary (§.= -1), compared to the eight node ele- 
ment HISQUE. Since addition of two nodes per element used on the hole 
surface does not increase the total global degrees of freedom significant- 
ly, no static condensation is performed by element HEL710. 
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5.3 Formulation of Element HEL710 


5.3.1 The Assumed Stresses 

With all the nodes present, element HEL710 possesses 20 degrees 
of freedom and therefore requires 17 betas in the assumed stresses to 
exclude spurious zero energy modes. For completeness, the cubic assumed 
stresses of element HISQUE, equations (4.4), with 18 betas were used. 
Since element HEL710 incorporates five nodes on one edge, it takes ad- 
vantage of a quartic displacement interpolation on this edge and hence 
for a better balance a complete quartic stress distribution, employing 
25 betas, was also tried. A complete quartic stress distribution is 
obtained by addition of the following terms to equations (4.4): 

<r„ = Yfa i 4 fa t fa fa f tyfa + yfa 

(Tyy s 6 / J/ + zL J^>%! + 4/- y ' Y>Z4 / X / (5.11) 

(fry = - 4^y An . 1 icy 1 fa - j fafa 1 - 

5.3.2 Interpolation Functions and the Strain Matrix 

The strain matrix is determined in a completely analogous manner 
to that of element HISQUE, with the only difference existing in the 
interpolation functions. The interpolation functions tr were constructed 
such that the displacement distributions on any boundary of the element 
were compatible with the number of nodes present on that boundary. 

Another requirement satisfied by the interpolation functions is: 

Jli - 1 
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The following interpolation functions were obtained for the isopara- 
metric element: 


^ - 5 (H) 1 (H )< 1 - OO-i )-4 4 4 


^ s i ('-VO* D - 4 4 4 - 4 <7 *>1 
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It is seen that by employing the ^ function, any combination of nodes 
numbered 8, 9 and 10, Figure 11, can be removed without disintegration 
of interelement compatibility. This characteristic will allow element 
HEL710 to be used in analysis with four and/or eight node elements. 

5.3.3 Imp lementation 

A 5 x 5 Gauss numerical integration was found to be the minimum 
order for an adequate integration of matrices H and <5. Very insignif- 
icant changes of the stiffness matrix resulted when higher order inte- 
gration was used. Otherwise, the implementation of element HEL710 is 
identical to that of element HISQUE, presented in Chapter 4. 

5.4 Results and Discussion 
5.4.1 Plate with Circular Hole 

The plate was modeled entirely with displacement based eight node 
elements ELEM8, except around the hole where elements HEL710 were used. 
Two sets of stress solutions were obtained. The first set of results 
acquired using complete cubic assumed stresses was completely unaccept- 
able. Stress concentrations of 1.80 a ) and 3.43^? were obtained at the 

* 0 O 

point of maximum stress, with Mesh-1 and Mesh-2, respectively. With a 
complete quartic stress distribution the stress solutions were somewhat 
improved, and stress concentrations of 3.03 and 4.21.0 were computed. 
The percent error in the stress solutions are summarized below: 
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Figure-12 Maximum Stress 3y HEL710 In 
Plate With Circular Hole 


(DOF) 




Element 

Stress 

% error 

Mesh-1 

Mesh-2 

Cubic 

-58.3 

-20.6 

Quartic 

-29*9 

-2.55 


Although not as accurate as element HISQUE, with Mesh*- 2 and a com- 
plete quartic stress distribution the computed stress concentration fac- 
tor converged within the +5% acceptable margin* The convergence of the 
stress solution is shown in Figure 12 for both cases. It is clear that 
element HISQUE with cubic assumed stresses is not effective in analysis 
and was therefore not used for the elliptical hole. 

5.4.2 Plate with Elliptical Hole 

Again element HEL710 was used around the hole while the rest of 

the plate was modeled with ELEM8 . Stress concentrations of 4.76^0 and 

8.47.40 were computed at the point of maximum stress with Mesh-1 and 

Mesh-2, respectively. The analysis did not converge to the correct 

solution of 9.50./? and proved element HEL710 unsuccessful, in all 
'0 

respects. Firstly, element HEL710 employs 25 betas resulting in a matrix 
of order (25 x 25) . Inversion of this matrix makes element HEL710 
almost twice as expensive as element HISQUE. Secondly, element HEL710 
employs four degrees of freedom more than element HISQUE, resulting in 
a larger global stiffness matrix and hence additional computational 
cost. 

5 . 5 Concluding Remarks 

The normal stresses around the boundary of the elliptical hole, as 
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determined by element HEL710 with Mesh-2, are plotted in Figure 13, 
against the polar coordinate location of the hole 9 . As angle Q in- 
creases, the ellipse exhibits more pronounced curvature variations; 
subsequently element HEL710 is unable to contain the magnitude of the 
normal stress close to zero. Elements inaccuracy is also partially due 
to the unbalanced nature of this element. Introduction of additional 
nodes will increase the cost of the element and further deteriorate 
any possible ability to compete against element HISQUE on economical 
grounds . 

Development of a more balanced element (i.e., comparable number of 
nodes on boundaries) which identically satisfies the traction free con- 
dition at the hole is judged most effective for analysis. 
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CHAPTER 6 


A Four-Node Element With a Circular Traction Free Edge 
6. 1 Introduction and Objective 

Instead of pointwise satisfaction of the free traction condition 
at the hole, via boundary point matching or other similar techniques, 
it is possible to formulate an element which will exactly satisfy this 
stress boundary condition at a circular hole. To this end a four- node 
element which incorporates a built-in traction free circular edge is 
formulated in polar coordinates. Although restricted to analysis of 
circular holes, the four-node polar element, named PET4, will serve as 
a model element for comparative evaluation. 

Element PET 4 is formulated with two different sets of assumed 
stresses. The first set of stresses satisfies the equilibrium conditions 
while the second set satisfies the stress compatibility conditions in 
addition to equilibrium. The objective is again for the improvement of 
solution accuracy, and efficiency, compared to other efficient alterna- 
tives, namely the mixed model approach of Chapter 4. 

6.2 Formulation of Element PET4 

The formulation of element PET4 is based on the modified complemen- 
tary energy functional / 7T MC . In the formulation the matrix Cj is 
determined by integrating the product of boundary tractions and inter- 
element displacements, around the element boundary (equation 2.7) . Since 
element PET4 has prescribed zero tractions on the circular edge, evalua- 
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tion of matrix <j is needed on the remaining three edges only. Figure 
14 shows the nodes of element PET4. 

6.2.1 The Assumed Stresses 

If the origin of the system of coordinates is located at the 
center of the circular hole, then the condition of zero traction at the 
hole reduces to simple constraints imposed on the stresses in polar 
coordinates : 


' rr 


ra 


\r=<z 


= 0 
r m a 


where: 

CL 5 Radius of the circular hole. 


(6.1) 


A set of assumed stresses that identically satisfy the equations of 
equilibrium can be obtained from the Airy stress function [14] . 


<Tr d 
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r dr r l ^e 2 - 

d 7 -# 

<?r z 
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Assuming the Airy stress function to be of the form 


( 6 . 2 ) 


d> - A r i (di'd f (J ox. l 

F /C n ) { 


(6.3) 
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Figure-l4 The Polar Traction-Free 
Element PET4 ~~ 
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then substitution of equation (6.3) into equations (6.2) yields: 
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Expanding. equations (6.4) with both terms and will provide 

a set of stresses weighted with respect to the coefficients • Im- 
posing the zero traction condition of equation (6.1) will reduce the 
number of constant coefficients and hence will in general provide an 
infinite series of polynomials in r and transcendental circular func- 
tions of 6 . An expansion of the terms with: 

1=0,1 

ft =0,1, Z, 3, 4, 5 


resulted in a twelve- term set of stresses. After rearranging the terms 
in increasing powers of r and redefining the coefficients * 

the following result is obtained: 
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or in matrix form: 
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PS 


( 6 . 6 ) 


Equations (6.5) show that retaining all 12 betas in the assumed 

stresses will result in a cubic distribution* while truncating the 

series af ter will render a quadratic distribution of the stresses. 

/ * 

The full set of stresses of equations (6.5), with twelve betas, was 
used in element PET4/12, and the truncated series with nine betas was 
used in element PET4/9. Equation (2.6) will then determine matrix H . 
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6.2.2 Displacement Interpolations 


The displacement interpolations are required along the boundary 
for evaluation of the £7 matrix. Along the edges connecting nodes 1 — 2, 
2-3, and 3-4 (Figure 14) , a linear interpolation between the nodal dis- 
placements will ensure interelement compatibility. 

Along any straight edge i connecting nodes i and i + 1 (Figure 15) 


Hi 



(6.7) 


It should be noted that displacements are interpolated in the Cartesian 
coordinates to finally render the stiffness matrix in the Cartesian 
rather than polar system of coordinates. Therefore it follows that at 
edge i : 


# - | ( a u,- + ? t *;) 

Z ~ 


( 6 . 8 ) 


where : 




l ' z 

/6 ~ Usrftfv Aide x f’rcm nca’jL x . 

Equation (6.8) applies to displacements U* and Uy both. 
^ is defined: 


If the matrix 
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L 

Then from equations (6.7), (6.8) and (6.9) it follows for side i: 


*"i,U - j (■ '- ^ ) 
l-i,li,2 s i C 1 * ’ z - /, 2, 3 


( 6 . 10 ) 


Again, no interpolation along the circular edge (4-1) is needed since 
no integration is required at this edge. 

The ^matrix is then determined from: 




L dS 


(2.7) 


since l* is the matrix of displacement interpolations in Cartesian 
coordinates, a transformation of the stresses from their original polar 
coordinates is necessary: 
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( 6 . 11 ) 
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" ( 6 . 11 ) 

X s ‘Trasisfcr/wfi/reh rPdh/X 

O 5 Po/o.r btXStJ Stresses - PA 
~ fz. 

hence : 

P - TP 

~ ^ ~ (6.12) 

Where Z" is the familiar transformation matrix due to rotation of the 
axes : 

W ^ rt 1 _ Zrfirt 

W 1 Zrtirt (6.13) 

/Wti - mrj /?7 1 -n 2 ' 

dV - Cos e 
ft - S/n & 

& = Ang/z, oj£ ro/zubon 

In equation {2.1), y is defined as the matrix of direction cosines of 
Figure 15 : 

y* 0 Yy ] 

(6.14) 
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which when multiplied by the stresses produce the traction vector. 
From (2.7) : 


6 



L dS 


From (6. 12) : 


6 : ( P r T T y L ds 

J 0Y 

6 * f p r y L dS 


(6.15) 


From equations (6.13) and (6.14) it follows: 


v s T r y 


V 2 


m z y y + m y y r? z Yy + Va 


„2y K - flW Vy 


m Z Vy - nw >> 


(6.16) 


Zmn y x + y y ( r n in 1 ) y y + Y y (,n z . v z ) 


As defined by equation (6.16) , the matrix y will transform stresses in 
polar coordinates to tractions in Cartesian coordinates, at any point 
Q on a boundary. Thus all the matrices of equation (6.15) are defined 
for integration. 
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6.2.3 Implementation and Numerical Integration 

For a Gauss-quadrature numerical integration of equation (2.6) a 
coordinate transformation is necessary: 


+• *! 


" * tj F(r,»)rdr<U * 


(6.17) 




where : 


F(r,t) s P r S P 
X - JcLcad/ast /vao''? ■h’arisfarmaJre/i 


For the transformation of equation (6.17) the polar coordinates are to 
be expressed in thejf,^ system: 


& 


" fa.q) 


(6.18) 


r 




(6.19) 


From Figure 14 it is seen that if sides connected by nodes 1-2 and 3-4 
are radial, and hence functions of r only, then function f of equation 
(6.18) constitutes a linear interpolation of angles labeled 6 and 
in Figure 14. 


e , i ( 


#3 -0. 


f \ 


(- 


'j. 0Z \ 


( 6 . 20 ) 


Where 5 is the non-dimensional coordinate of equation (6.17) defined to 
assume values of +1 and -1 at the boundaries 1-2 and 3-4, of Figure 14, 
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respectively. This assumption forces element PET 4 to possess two 
radial edges, but its adoption results in considerable simplifications 
and is hence preferable. 

With $ interpolation fixed, interpolation in the r direction is 
straightforward. Defining the non-dimensional coordinate // to be -1 at 
the circular edge 4-1 and +1 at the straight edge 2-3, then: 

r = a i /j S - 1 ( 6 . 21 a) 


r r h(4) \ = 


(6.21b) 


tihlrti Ml-tjYjL ( 6 . 22 ) 

fXz-X . - (iz-Wb** 

is the equation of line 2-3 in polar coordinates. In equation (6.22), 

x , y are the Cartesian coordinates of node i. The interpolation in 
i i 

the r direction is then simply: 


r s i [Jue) - 


(6.23) 


It is seen that equation (6.23) reduces to the conditions of equations 

( 6 . 21 ) . 

With functions f and g of equations (6.18) and (6.19) defined the 
Jacobian of the transformation is given: 
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From equation (6.20), it is clear that 9 is independent of and there- 
fore the determinant of the Jacobian matrix simplifies: 


+ fu*) | (6.24) 

For a numerical integration of the matrix/^ as shown in equation (6.17), 
Gauss integration stations ? ^ and ^ are picked and their corresponding 
r and 9 coordinates are determined from equations (6.20) and (6.23) . 
This in turn enables the evaluation of the stress matrix P at any inte- 
gration point. Equation (6.17) then reduces to the familiar form: 


H t t 


Z 


j 

y 

♦ 

/ 


F (i'tfj ) 7 (l"Tj )j W* sJj 


(6.25) 


where : 

= Gauss integration weights. 

For the numerical integration of equation (6.15), the linear interpola- 
tion of equation (6.8) is used: 



X C^ft ) « - X* ) + ~ (*’*■-! + } 


yd*) s ~r (%h ~y*) + t 


; £<?££ (i,ifl) 


where: £ B 6sluS$ /rthfiraXcn iX&han 

*) 0 


therefore : 


r(C„) ; YXWi 

’ya«) | 


Or 0 (in ) - 


/O-rt 


/■<*>. 


Direction cosines of equation (6.16) are also defined for edge 

(i, i + 1) 


Ys 


X " 


V, 


r * 


Cos 1 

I ~ X- 4 - 

l \ y - , Vi H 

Sjo | 

r T'i a ( X 4 *■ 

t 

1 /,• - /-V/ 


Equation (6.15) can then be broken into three separate integrals along 
the three straight edges of the element 


/S _ 


(6.26) 


A S- 
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where : 



(6.27) 


From equation (6.15): 



And from equation (6.8) : 

da _ / 

di I 

Equation (2.18), along with (6.25), (6.26) and (6.27) determine the 
stiffness matrix j(. . 

In most cases a 5 point Gauss integration was found sufficient for 
accurate integration of matrices// and . For the complete stress 
distribution of equation (6.5) with 12 betas, 7 points were required to 
integrate the 0 matrix. 

6.2.4 Element PET4X5 

As previously discussed in Chapter 4, in two dimensional plane 
analysis, the assumed stresses can be chosen such that the condition 
of compatibility is also satisfied: 

y '*> - o (4.i5) 
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If the Airy stress function is assumed to be [14] [15] : 

oo 

f = f 6> <r) f 2L C r ) i^/d/unr>&‘, d/JxifiS \ 

= 4,-4 hr , C . r'-f d . r z J.„r 


J (/)) . /> / _-/7 Z+ ft J 

(r ) s & v F /. bn f + C# r d- a, ? r 


2-/7 


(6.28) 


i ( 0 ) 4 . 


then substitution of ^ in equation (6.2) and (6.1) gives: 

cd • b > ( j 2 - /J *■ ld ° ( tnr - l " a ) 

Tfia - — Sa ( — ^ - 4^ j f 2 do / />; r - hi A ^ 


rr (0) 

&r4 ' 0 


(6.29) 


( n) 

Similarly by substitution of ^ it follows: 

(n) i 2 fi'Z -2n+2- fl-L 

V tr \ \(xK,)C n ar r 

v ln+Z ~/i'l 2 , v *2- -v-l 

-(/+n)C*<X r ~(/-n L )d,7d r 

f (Z+n-n z )Cn (Z-n- n z )dv r ' n daoas\ 


(n) ( 

04 ~ \ 


i _n-i 


-ZJt+Z /i -2 


C TS' * ) (/-n z )C«a~r n f (t-n)dn & r 


+ (/+n) Cn A 


ZfljZ - Z- n . 7. j . ~ A 

r + €7/ a f 


2 „ -*-Z 


- yy \ 

/ (Z r dn+/i z ) d/: a / + /Z. 2 si* n y d/?r ' fjC.h r .a . 


(6.30) 


/I 


■ ' 1 * •-■ 'Jt 

F0i& *V*iv 
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Vf*' -- j (l-n'^Cn&'r"' 1 * (l.n)dn a~ ln *~ r "' Z 


2 / 1+1 


- (/■ftl)c„ a -" r ~r (/. r'”' 2 - 

-r ( /, 'i+niCn /~ n +(n-n i ) dn r' n \ (. C*°n9 ; fiv*.ns\ 


(6.30) 


In condensed matrix notation: 


CO 


<r c t ( °\ Z * 


Ci ) 


n 


.(0) 


If the logarithmic terms appearing in £T are dropped ; for n 
the following results are obtained: 

Vrr - 0 -~z )/-' t( r 'fj) { t/L*e/U j 

f (If i£‘ l ) l C/Mufr tsd^wfs) 

r it , f(Zr t £?) 

_ ( / j . zZL \ ) ' 2 & A — xA ’/'^ s / 

‘ /~4 1 | / ' j 

tf r0 r J - Cc^G'il \ 


= 1 , 2 


/'J / ! 


/ z _ >ki'_ / ) ' 

( - a- ' 


^ ~ ^ ( 

* J 

/i * c / „ /?. * ^ r* 


/. 


(6.31) 


where the coefficients C and D of equations (6.29) and (6.30) are 

n n 

replaced with „ ’ s. 

/ 
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Using the assumed stresses of equation (6.31), element PET4 will satisfy 
conditions of stress equilibrium, compatibility and free traction at the 
boundary, in addition to interelement compatibility. As shown in equa- 
tion (6.31) only 5 betas are present in the assumed stresses and there- 
fore the size of the matrix H is kept to a minimum. The element using 
the set of stresses of equation (6.31) is named PET4X5 and with the 
exception of the assumed stresses, its formulation is identical to 
elements PET4/9 and PET4/12 discussed earlier. 

6 . 3 Results and Discussion 

On the finite element model of the plate with circular hole ele- 
ments PET4 were used at the hole boundary. Since elements PET4 possess 
two nodes at each edge, a compatible four-node displacement based iso- 
parametric element [12] , called ELEM4, was used to model the plate 
away from the hole. 

The plate was analyzed with elements PET4/9 and PET4/12 to deter- 
mine the optimum set of assumed stresses between the two. Element 
PET4/12 was found superior in performance although a converged stress 
solution was not obtained with Mesh- 2 in either case. Figure 16 shows 
a plot of the stresses obtained with the two elements. The percent 
error of the stresses are summarized in the following table: 


Element 

% error 

Type i 

VIesh-1 'Mesh-2 \ 

PET4/9 

| -30.7 i - 12 . 2 ! 

PET4/12! - 36.0 ; - 8 . 31 : 


OF F 


OUALry 
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Although the results may at first appear unsatisfactory, a closer 
comparison of Figure 16 and Figure 6 will indicate that analysis with 
elements PET4 utilizes less than half of the number of degrees of free- 
dom of an equivalent model using eight node elements. Element PET4/12 
is in addition much more economical than element HISQUE, when employing 
8 nodes and 18 betas. Combined with the four-node displacement based 
elements ELEM4, a reduction of more than 50% in the total solution time 
(CPU) resulted in comparison with the eight node analysis of the same 
mesh. 

Element PET4X5 is yet more economical than element PET4/12 in that 
it only uses 5 betas. With Mesh— 1, employing a total of 16 degrees of 
freedom, PET4X5 produced a stress concentration of 4.52 <p , already 
within the acceptable margin of ^+5%. With Mesh— 2 a stress concentration 
of 4.13 (error = -4.4%) was predicted. These results are also plotted 
in Figure 16. The in plane stress along the y axis, obtained with 

PET4X5, is plotted in Figure 17 for Mesh-1 and Mesh- 2. Element PET4X5 
is judged so far the most efficient and effective finite element in 
analysis of stress concentration around circular holes. 

6.4 Concluding Remarks 

Satisfactory stress solutions of element PET4X5 suggest that in 
addition to assuring the traction free boundary condition at the hole, 
satisfaction of stress compatibility conditions will further improve 
the performance of the element. In the case of element PET4/12, where 
compatibility conditions are not met, the stress solution contained 
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8.3% error when modeled with Mesh-2. In view of the considerable cost 
reductions achieved with element PET4/12, in addition to the fact that 
few degrees of freedom (42) were used in analysis, this element presents 
desirable advantages compared to the eight node hybrid isoparametric 
element HISQUE. The apparent failure of element PET4/12 is nevertheless 
more due to the ineffectiveness of the four-node elements in connects 
than its own capabilities. 
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CHAPTER 7 


Conclusion 


7.1 An Eclectic Approach 

Use of hybrid isoparametric finite elements was found more desir- 
able than their displacement based counterparts in stress concentration 
analysis around holes. Although these elements do not strictly satisfy 
the required traction free boundary condition, they consistently approx- 
imated this condition more accurately than the conventional displacement 
element. 

Adoption of the boundary point matching technique was seen to reduce 
element cost effectiveness without a justifiable proportional increase 
in element accuracy. However, improvement in solution accuracy and 
efficiency was achieved when the free traction condition at the hole 
was exactly satisfied by adoption of an appropriate set of assumed 
stresses. Here, hybrid formulation provides the only effective means 
of satisfying stress boundary conditions. 

Further improvements resulted when in addition to satisfaction of 
boundary traction, stress compatibility conditions were also satisfied. 
This deduction is a natural consequence of satisfying stress equilibrium, 
compatibility, boundary traction and interelement displacement continuity 
in such finite elements. 

In three dimensional elements, satisfaction of stress compatibility 
may prove too difficult although other requirements met in two dimen- 
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sional analysis can be assured. A general formulation of two dimension- 
al hybrid isoparametric elements that incorporate any number of traction 
free boundaries is provided in Appendix A. This formulation can be ex- 
tended to solid elements and its generality affords satisfaction of 
free traction condition on boundaries of arbitrary geometry . 

Finally, employment of such specialized hybrid elements leads to 
an eclectic modeling of the physical problem which in general yields 
more efficient and accurate solutions. This implies use of specialized 
elements where stress boundary conditions are prescribed, hybrid elements 
in regions of high stress gradients, and displacement elements other- 
wise. Indeed in problems of stress concentration, fracture mechanics, 
or wave propagation use of hybrid elements will provide noticeable im- 
provements . 

7.2 Recommended Research 

The twelve node solid element of Figure 18, as well as the two 
dimensional six node elements it incorporates on each, side can be for- 
mulated by the approach of Appendix A. Performance of the six node 
element should be compared with elements HISQUE and PET4 of Chapters 4 
and 6. In this element node 5 is introduced on the mid-boundary to 
provide a traction free edge of flexible geometry, while node 6 is to 
allow connection with eight node elements for added effectiveness. 
Satisfactory results with the two dimensional element should justify 
its extension to a study of the twelve- node solid element. 
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APPENDIX A 


Satisfaction of Free Traction Boundary Condition in 
Assumed Stress Hybrid Isoparametric Finite Elements 

A. 1 Introduction and Objective 

In finite element analysis of stress concentrations induced by the 
presnece of holes in a structure, satisfaction of the free traction 
boundary condition at the hole becomes a prerequisite for effective 
stress computation. Since isoparametric finite elements conveniently 
conform to shapes of arbitrary geometry, formulation of a two dimension- 
al hybrid isoparametric element that identically satisfies the condition 
of free traction at the boundary is developed. The general characteris- 
tic of this formulation makes it applicable to three dimensional solid 
elements as well. 

A. 2 Curvilinear Coordinates of Isoparametric Transformation 

In isoparametric transformation, Cartesian coordinates X and Y 
are interpolated by a set of shape function hr : 







y* 



(A-2) 
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QUALITY 



where: 


(X., Y.) » Coordinate of node i. 

1 1 

fl • Number of nodes. 

Equations (A-l) and (A-2) can be rewritten in summation form: 



(A- 3) 


(A-4) 


In equations (A-l) through (A-4) , h^ are functions of the isoparametric 

_ / t 

coordinates £ and g . From these equations the covariant base vectors 
of the isoparametric system are found [16] : 

l ■ ( *t f t , ) ** ( y f % ) ** 

& • (*t %Y'* (' J f %)* 

from which the covariant metric tensor results: 



(A-7) 
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(A-8) 


L - <*f 

6 - **<&) 



(A-10) 


The normal vector to the plane of the element, /Ip can be determined 
from: 

n f - . is 

f /l 


The contravariant base vectors of a general isoparametric planar ele 
ment follow: 



The Christoffel symbols of the transformation are defined: 




O £ 

V pm 

2l n 
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Defining: 


3/ 

■dv 



3Y 

H * 








C - -i 

ft 












(A-13) 


(A-14) 


(A-15) 


( A-16) 


(A-17) 


(A-18) 


A. 3 The Boundary Normal and Traction Constraints 

Consider the eight node isoparametric element shown in Figure A-l, 

with vector T designated as tangent to boundary ^ = -1 . This boundary 

constitutes a parabola, the equation of which is determined by element 

shape functions h^, h . and h . 

2 3 6 
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r 3 (A 2 ^2 + ^3 ^3 -t Xt ) 'Ll + ( ^2 ^2 v Si ^3 + Vi, r^6 ) 


(A-19) 


where : 


Hi 


k 



Hi (t z ) 


The tangent vector T follows directly from (A-19) : 





'dHi 

n 2 



(A-20) 


If the tangent vector is defined: 

7 ” = ( A -Lj + S 'i.z ) 

then the normal vector to this boundary is: 

N =s (S 'Ll - A iz ) (A-21) 

N ' T * 0 

It should be noted that 7" , and /V are not unit vectors, but can be nor- 
malized if needed. From equations (A-20) and (A-21) the normal vector 
results : 


88 



(A-22) 



Equation (A-22) is then compared with the first contravariant base 
hector given by equation (A-ll) . At the boundary Jr .1 : 



Noting that the shape functions h^ are defined to vanish at nodes i / g, 
equation (A- 23) reduces to: 



The vectors Si and ^ defined by equations (A-22) and (A-24) are seen 
to be parallel since: 





Therefore the normal vector at the boundary 5 ^ is parallel to the 

first contravariant base vector A . As a result of this conclusion it 
can be shown that if the assumed stresses used in the hybrid formulation 
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are the contravariant tensor stresses of the isoparametric system, the 
condition of free boundary traction reduces to simple conditions of con- 
straint on the stresses. Similar conclusion can be drawn by considering 
the boundary ^ L +1 . 

Along the boundary £. = constant, components of the traction vector T 
are expressed in terms of the contravariant stress tensor: 


t 



r 


y 


n 


J 



r ,i n i 


L /?. 


(A-25) 


(A- 26) 


f 




(A- 27) 


Where i •' is the contravariant stress tensor, and n i are the components 

A 

of the unit vector N, normal to the boundary: 


A 


/V r 





(A-28) 


The normal vector N was shown to be parallel to the contravariant base 
vector ^ and therefore: 




(A-29) 
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Comparison of equations (A-28) and (A-29) gives: 


n, - JL 
/g" 

rii * o 

which reduces equations (A-26) and (A-27) : 

t ' , T"n, 

t 1 T 11 '!, ^ t ll n, 


The traction vector is resolved into contravariant components 
covariant base vectors: 


t , 



+ 



substituting from equations (A- 30) and (A-31) : 

T * T N 'll + £ 11 n / ^ ; ("I. 1 = CanSnU^f ) 

, / 

For zero traction at the boundary g. = constant: 

f r". r' 1 - 0 ] , 

[ J £ — CO r\£-trxn+ 

Similarly for a prescribed pressure d q at this boundary ; 


(A-30) 


(A-31) 


along the 


(A-32) 


( A- 3 3 ) 


(A- 34) 
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which is seen to reduce to equations (A-34) in the special case of no 
applied pressure. 

* 

Similar conditions are obtained at boundaries constant. As 

shown in Figure A-2, the first covariant base vector ^ is always 

tangent to this boundary. From equation (A-12) it follows that the 

2 2 

contravariant base vector £ is normal to line i % constant. The unit 

D ^ 

normal is hence obtained: 


/< 


A/ r 







(A-36) 


Comparison with equation (A- 28) shows that along lines 



constant: 


n, eO 


(A-37) 


Substituting (A-37) into (A-26) , and (A-27) will in turn simplify equa- 
tion (A- 32) : 






\ 


For 




at boundary > - 


constant : 


OH'-mil page is 
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0 


(A-38) 



as t 


1 


A. 4 Implementation 

The following formulation illustrates how the aforementioned re- 
sults can be applied to a two dimensional isoparametric element incor- 
porating a traction free boundary at r . 

A. 4.1 The Assumed Stresses 

The two dimensional contravariant stress tensor can be obtained 
from the Airy stress function: 






(A- 39) 


.f 


where £ is the two dimensional permutation tensor: 


e 


// 



as 0 


ft 


and £ is defined by equation (A-10) . 

O 

Stresses obtained from euqation (4-39) will necessarily satisfy 
the homogeneous equations of stress equilibrium: 



- /) 
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From equation (A-39) it follows: 


r" * 

/ ^ ?-<S 

5 {(n 1 )r ~ W 

r? 7 

u 

</ r 

r* i 
" J 

(A-40) 


/ l p-V 3$ 

5 l?/';- ' ?i’ 

P ! _ 

u 

£ ‘r* 

P<rZ 

r zi , 
" j 

(A-41) 

v 1 . 

_ £ j 7£ r 7 _ ££ 

£ 1^7 ^ 

; 

Vi 1 

J 

(A-42) 


; rjfl s ChrishfftJ Symbols 

For simplicity the Christoff el symbols are redefined: 


r u 

&a,n) = 7 ; 

M ~! 

tZ 

a 

* r 1 

/i 



D(l'H) 
£ 7 /■: 




// 


With. ? and r' replacing J and g "respectively. Substituting into 
equations (A-40) , (A-42) , and (A-41) gives: 
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f "* T l f,i i - 4>,1 \ 

b 1 ' ; 

I . £ U 1) 4>„ - 0(1, n) £. | 


<5 


T lt =, Z | <*,, - £G,p 4„ - Hi'ftbi ] 

t 

where C and are partial derivatives of <fi with respect to | and tj 
respectively. 

Assuming the Airy stress function: 


d> = X7i Ann l n J 

l /7i n 1 


the contravariant stress components become: 


? . i it < (/i-t Z) (n+i ) Am n+z 
8 * 


(A-43) 


/? 

r . 


- < ("to A„ fh „ Aa,f ) - &(i‘i ) j * 

-(W-t 1 ) A/r]^/ t n C > ,// / ) 


(A-44) 


(>ir ! ) A™ , *t' Dili*) \ 
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- (™t> ) 

1 (A-45) 

- (fl+‘) Am, /it! A(lt^) j 


To prescribe zero tractions at boundary 2 s .| of the element: 


r" = t /z . o 




(A-34) 


Substituting equations (A-43) , and (A-44) into equation (A-34) introduces 

constraints among coefficients A . Rearranging the indices to unify the 

mn 

coefficients A gives: 
mn 


M A/ 


M A/ 


T a Z,JL Ini*"*/)* C t * i'**!*' 1 5 \ A™* 

£ /».***•{ l c / j 


( A- 3 5 ) 


T“ , I ZZ ft l ') n ' 1 

£ mzc n~o C L i j 

1 ZZ { *(«-» "£ a,?) (a-36) 

- n i * y*' 1 F(l,y) | A™ 


(A— 37) 


where the functions A, B, C, and 5 are defined: 


6(i> 1) 
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and therefore vanish at . As a consequence T and T of equations 

(A-35) and (A-37) are also identically zero at this boundary. Equations 

(A-35) through (A-37) can be expanded in as many terms as required. 

Redefining coefficients A as A; , the stresses are expressable in 

mn J * 

matrix form: 


£Y i >]) ~ /l 


(A-38a) 


An advantage of this approach is immediately obvious from equation 
(A-38a) . Since the assumed stresses are chosen in the isoparametric 
system, the resulting stiffness matrix is invariant and hence unaffected 
by element rotation relative to the Cartesian coordinates. This elim- 
inates the need for a complete polynomial of stresses and allows selec- 
tive adoption of effective terms. 

For implementation on the high speed computer equations (A- 38) can 
be rewritten in a much simpler form: 

O' = r p /■ (A-38b) 

which compared with equation (A- 38a) gives: 





r p 


(A-38c) 
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Here V is the matrix of Christoff el symbols, the entires of which were 


previously defined: 

0 A £ 

/ - £■ -P (A-38d) 

0 , C -0 

Matrix P of equation (A-38c) is simply obtained by expansion of the 
following series: 


r s i 

0 


M A/ 

/ . zz 

/><zo rtze 


m-l ti 

*r i"'' 


(A-38e) 


Thus equations (A-38c) , (A-38d) , and (A-38e) completely define a set of 

self equilibrating assumed stresses. 

In computation of matrix fZ : 


0 r f P r s P a'V 


( A-39) 


Gauss integration coordinates become the natural variables and the neces 
sity for evaluation of a Jacobian matrix of transformation is eliminated 
The assumed stresses of equation (A- 38c) can be easily evaluated at 
Gaussian integration stations and equation (A-39) becomes: 


//-- 1 ZZ ( p t sp ).. 


“ti 


u / * u: # 

r/ / ' • / 


///, n ; \ ? ■ZAo'f .v.rt/rpr’/: 


■* Ui* 

• .v » 


l - ' 



( A-40) 
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The compliance matrix 5 of equations (A- 39) and (A-40) relates the co- 
variant strain tensor to the contravariant stress tensor: 



( A-41) 


For isotropic materials 5 is defined [17] : 



>5 s Young's modu/us. 
y s Potsson's tfaha . 


(A-42) 


In matrix form equation (A-41) becomes: 


J// 


5//// 

S//zi 

l S///z 


r" 

£,22 

> S 

Szju/ 

Szzzz 

jLSzz/l 

< 

r u ► 

i/z 


S/zu 

S/ZZ2 

zsm 


r /2 
► « 


wher ' : S //// = i 

£ On 

5// “ • ■ j { < ur >6* - v 6„ 6 * ) 

S/ "i * Swi ' j6„ 6,i 

S Zl U * 1 £ 

Szza -- j 6zz • S a12 - 

5 -iJ { 6zz * ] 
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A ,4. 2 The Strain Matrix 


Interpolation of the strain matrix 3 is required for computation 
of matrix 3 . In the Reissner formulation: 


6 , 

S'* 

To obtain the final stiffness matrix in the Cartesian system of coor- 
dinates, displacements are interpolated in Cartesian coordinates and 
hence transformation of the stresses to the same system becomes neces- 
sary: 


(ft 


P r 6 dV 


(A-44) 


where : 


(f s vjf t n 

*t g con/ravasi<v>f fSapasanu/nc shisscs 



etudes ia-n 5 fosses (XiV) 


In matrix form, equation (A-45) becomes: 


1 


’ Tn 

T* 

7/3 


'T"' 


► » 

Tv 

Tiz. 

7z3 

4 

T u ► 

<r>i 

m « 


T3, 

7~3Z 

T33 


r l 


(A- 46) 
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where : 


r„ - 


T* 


(M) z 


• (£)' 


r, 3 * 

Tu * 


n 

2 . 2* 21 

31 3i 

<&• 


Tzl = (21 \ z 
v- dn / 

TEs - z 2/ 

?j ?ri 

75/ * 2A M 

P£ 7*1 
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31 *1 

T 33 •z2&21 

31 If 


Derivatives of the Cartesian coordinates with respect to £ and 
are determined from the interpolations: 


4 

y * tji fa 

jk. 


(A-l) 


(A-2) 


In isoparametric elements the displacements are interpolated identically 


IL Zj U A hi 

X 

IT* ZflTihi 


(A-47) 


(A-48) 
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From equations (A-47) and (A-48) the strains can be determined: 
<£ // r 




= 2jl 

W 


r 2 ^ . %£ 

?y i ' dt- 

These equations in turn determine the matrix & in terms of the nodal 


displacements 





(A-49) 


With equations (A-46) and (A-49), matrix G of equation (A-44) becomes: 



T P 5 dV 

/v 


(A- 50) 


where "7" is the matrix of transformation, and P is defined by equation 
(A-38c) . Matrix G of equation (A- 50) can also be integrated by the 
Gauss quadrature scheme of equation (A-40) . 


A. 5 Extension to Solid Elements 

Similar results can be derived for solid isoparametric elements, 
Figure A-3. The three dimensional components of a traction vector 
expressed in terms of the contravariant stresses are: 


t 1 - t" n, t ? a r)i t f' 3 />3 

t* » t v n, * , T i3 n 3 
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J. 3 + 1‘ + H *33 „ 

t * T r h + T ni + T ^3 


(A-51) 


where : 


A 

» Components of- unit fidroiaJ tftcJor A/ 


At the boundary £/•-/ of Figure A- 3, the contravariant base vector 
^ is along the direction of the normal to this boundary, A/ : 


A 

N 


A/ 

/A// 



; fll s fli s O 


(A-52) 


Substitution of equation (A-52) into equations (A-51) gives 

t ‘ « fn. 

t 2 r T 1 ' ni 

t 3 = t 1 ' m 

Therefore if the assumed contravariant stresses of a hybrid solid element 
are chosen such that: 



the condition of free traction will be automatically satisfied at bound- 
ary I — of this element. 
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APPENDIX B 


Computer Programs 

Ttie conputer programs are all in form of subroutines. This 
enables use of the programs in Finite Element codes by way of 
linking a main program to the desired element stiffness matrix 
subroutine. Ehch program is followed by another subroutine for 
computation of the stresses. 

All subroutines were programmed and executed on the Digital 
VZQC-1 1/780 conputer. All programs are written in FORTRAN-IV. 


B.l ELEMENT HISQUE 

The following variables are required by subroutine HISQUE: 
E =Young 1 s modulus 

NU =Poisson’ s ratio 

T =Element thickness 

X =Vector of nodal coordinates, x 

Y ^Vector of nodal coordinates, y 

! x =0(1) if node x is absent (present) 

NBETA =Nunber of betas 
10 Output control variable 

The following variables are returned by subroutine HISQUE : 
HINVG =Product of matrices [H]inv. and [G] 

REALK =Element stiffness matrix 
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Subroutine QUEST conputes the stress vector STRESS f rcm the 
following variables: 

XS, YS =Coordinates of point at which stress is desired 
Hinvg =Matrix returned by subroutine HISQUE 
DISP =Vector of nodal displacements 
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c************** *********************************************** ******* 
KUROSH KAFIE MARCH 1981, M.I.T. * 

******************************************************************** 


SUBROUTINE HISQUE (E,NU,T,X,Y,I5,I6,I7,I8, NBETA, 10, HINVG, REALK ) 


IMPLICIT REAL*8 ( A-H , 0-Z ) 

DOUBLE PRECISION NU,J11,J12,J21,J22 
REAL HINVG(NBETA,1) ,REALK(2* (4+I5+I6+I7+I8) ,1) 
CHARACTER*5 ASTRNG 
CHARACTER*1 3 BSTRNG 

DIMENSION X(l) ,Y(1) ,GS (4, 4) ,GW(4, 4) 

DIMENSION DIM (3) ,F(8) ,DZ(8) ,DE(8) 

DIMENSION S(3,3),P(3,18),H(18,18) 

DIMENSION BE(3, 16) ,G(18, 16) 

DIMENSION HG (18,16), STMTX(16, 16) 

DDX(M)=-hJ22*DZ (M)-J12*DE(M) 
DDY(M)=^J21*DZ(M)-KJ11*DE(M) 

DATA GS/ O.OODO, O.OODO, 

1 . 00000000000000 ODO, .OOOOOOOOOOOOOOODO, 

2 -.5 773 502691 89626D0, .5 773502691 89626D0, 

2 .OOOOOOOOOOOOOOODO, .OOOOOOOOOOOOOOODO, 

3 -. 7 745 966692 41 483 DO, .OOOOOOOOOOOOOOODO, 

3 .7 745966692 41 483D0, .OOOOOOOOOOOOOOODO, 

4 -. 861136311594053D0, -. 339981 04358485 ®0, 

4 . 33998104358485 6D0, .861136311594O53D0/ 

DATA GW/ 2.00D0, O.OODO, 

1 .OOOOOOOOOOOOOOODO, .OOOOOOOOOOOOOOODO, 

2 1 . OOOOOOOOOOOOOOODO , 1 . OOOOOOOOOOOOOOODO , 

2 .OOOOOOOOOOOOOOODO, .OOOOOOOOOOOOOOODO, 

3 .555555555555556DO, . 888888888888889D0, 

3 . 55555555555555 6D0, .OOOOOOOOOOOOOOODO, 

4 .3 478548451 3 7454D0, .652145154862546D0, 

4 . 6521451 54862 546D0 , .34785484513 745 4D0/ 

NNDDES=4+1 5+1 6 +1 7+18 
NDOF =2*NN0DES 

C5 = 15 

C6 = 16 

C7 = 17 

C8 = 18 
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NINT=3 

IF(NN0DES.GT.4) NINT=4 
IF(NBETA.GT.12) NINT=4 

IF(NNDDES.LT.4.0R.NNDDES.GT.8)G0 TO 1000 
IF( (NDOF-3) .GT.NBETA) GO TO 2000 

IF(NBETA.GT.18) GO TO 3000 

DO 100 IR=1, 18 
DO 50 ICH=1 , 1 8 
50 H(IR,ICH)=0.D00 
DO 75 1 03=1,16 
75 G(IR,IOG)=O.DOO 
100 CONTINUE 

DO 145 IRGW=1,3 
DO 125 IC0L=1 , 1 8 
125 P(IROW,ICOL)O.ODOO 
DO 135 JCOL=l , 1 6 
135 BE(IROW,JCOL)=O.ODOO 
145 CONTINUE 

N5=I5* (5 ) 

NS=I6* (5+15 ) 

N7=I7*(5+I5+I6 ) 

N8=I8* (5+I5+I6+I7) 

IF (N5.EQ.0) N5=8 
IF (N6.EQ.0) N6=8 
IF (N7.EQ.0) N7=8 
IF (N8.EQ.0) N8=8 

S (1 , 1 ) =1 . 0D00/ E 

S (1, 2)=-NU/E 

S(l, 3)=O.ODO 

S(2,1)=-NU/E 

S (2 , 2 )=1 . 0D00/ E 

S (2, 3)=0. 0D0 

S(3,1)=0.0D0 

S ( 3 , 2 ) =0 . ODO 

S (3 , 3 )=2 . DO* (1 . D04NU) /E 

200 DO 500 II=1,NINT 
ZEX3S(II,NINT) 

DO 500 JJ=1,NINT 
ET=GS (JJ, NINT) 

OPZ =1 .DO+ZE 
OMZ =1 . DO-ZE 
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OPE =1 .DO+ET 
OME =1 . DO-ET 
CMZSQ=1 . DO-ZE* *2 
OMESQ=1.DO-ET**2 

F ( N5 ) =C5 *CMZSQ*OPE/ 2 . DO 
F(N6) =C6*OMESQ*CMZ/2 . DO 
F ( N7 ) =C7*CMZSQ*CME/2 . DO 
F ( N8 ) =C8*OMESQ*OPZ/2 . DO 

F(1)=0.25DO*OPZ*OPE-(F(N5)+F(N8) )/2. 
F(2)=0.25D0*CMZ*OPE-(F(N5)+F(N6) )/2. 
F ( 3 ) =0 . 2 5D0 * CMZ *CME- ( F ( N6 ) +F ( N7 ) ) / 2 . 
F(4)=O.25D0*OPZ*OME-(F(N7)+F(N8) )/2. 

XSO.ODOO 
YS=0 . ODOO 

DO 250 K=1 , NNODES 
X5=XS+X(K)*F(K) 

YS=YS+Y(K)*F(K) 

250 CONTINUE 

XSO=)S*5S 

YSQ=YS*YS 

XCU=JS*JSQ 

YCU=YS*YSQ 

P(l/ 1)=1.0D00 

P(1,4)=¥S 

P ( 1 , 6 ) =XS 

P(1,8)=YSQ 

P(1,10)=)SQ 

P(1,11)=2.D0*XS*YS 

P ( 1 , 1 3 ) =YCU 

P(1,15)=XCU 

P(1,16)=XS*YSQ 

P(1,17)=X5Q*YS 

P ( 2 , 2 ) =1 . ODOO 
P(2, 5)=XS 
P(2, 7)=YS 
P ( 2 , 9 ) =XSQ 
P(2,10)='iSQ 
P(2,12)=*>(1,11) 

P(2,14)=XCU 
P(2,15)=3.D0*X5*YSQ 
P(2,17)=YOJ/3.DO 
P(2, 18)=XSQ*YS 
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P(3,3)=1.0D00 
P(3,6)=-YS 
P(3, 7)=-XS 
P(3,10)=-P(l,ll) 

P(3,11)=-YSQ 

P(3, 12)=-XSQ 

P(3,15)=-3.D0*XSQ*YS 

P(3,16)=-YCU/3.D0 

P(3,17)=->5*YSQ 

P(3,18)=-3OJ/3.D0 

DZ(N5)=-C5*ZE*OPE 
DE(N5 )=4C5*OMZ9Q/2 . DO 
DZ(N6)=-C6*CMESQ/2 .DO 
DE(N6)=-C6*OMZ*ET 
DZ(N7)=-C7*ZE*CME 
DE( N7 ) =-C 7*OMZSQ/2 . DO 
DZ (N8)=4C8*CMESQ/2 . DO 
DE( N8) =-C8*OPZ*ET 

DZ ( 1 ) =0 . 2 5DO*OPE- ( DZ (N5 ) +DZ ( N8) ) /2. 
DE(1)=0.25DO*OPZ-(DE(N5)+DE(N8) )/2. 
DZ ( 2 ) =- . 2 5D0 *0 P E- ( DZ ( N5 ) +DZ ( N6 ) ) / 2 . 
DE(2 ) =0 . 2 5D0*0MZ- ( DE (N5 ) +DE (N6 ) ) /2 . 
DZ ( 3 ) =- . 2 5D0*CME- ( DZ ( N6 ) +DZ ( N7) ) /2 . 
DE(3 ) =- . 2 5DO*OMZ- ( DE (N6 ) +DE (N7 ) ) /2 . 
DZ ( 4 ) =0 . 2 5D0*CME- ( DZ ( N7 ) +DZ ( N8) ) /2 . 
DE(4)=-.25DO*OPZ-(DE(N7)+DE(N8) )/2. 

J11=0.0D00 

J12=0.0D00 

J21=O.ODOO 

J22=O.ODOO 

DO 350 K=1 , NNODES 
J11=DZ (K) *X(K)-KJ11 
J12=DZ(K) *Y(K)-KJ12 
J21=DE(K)*X(K)+J21 
J22=DE(K)*Y(K)-KJ22 
350 CONTINUE 

DET J=CT1 1 ’•U22 -J1 2 AJ 2 1 

WTH=GW (II, NINT ) *GW ( J J, NINT ) *T*DETJ 

DO 400 J=1,NBETA 
DO 370 K=l, 3 
DLM(K)=0. ODOO 
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DO 370 L=l, 3 

370 DUM(K)=€>UM(K)-«(K,L)*P(L,J) 
DO 390 I=J, NBETA 
HTEMP=0 . 0D00 
DO 380 L=l, 3 

380 HTEMP=HTEMP+P(L f I)*DLM(L) 
390 H(I,J)=H(I,J)+HTEMP*WTH 
400 CONTINUE 

DO 435 1=1, (NDOF-1) ,2 
BE(l,I)=ODX( (l+l)/2) 
BE(3,I)=DDY( (l+l)/2) 

435 CONTINUE 


DO 445 I=2,NDOF,2 
BE(2,l)=ODY(l/2) 
BE(3,l)=ODX(l/2) 
445 CONTINUE 


WTG&ti ( 1 1 , NINT ) *GW ( J J, NINT ) *T 
DO 450 1=1, NBETA 
DO 450 J=1 , NDOF 
DO 450 K=l, 3 

450 G(l, J)=<j(I, J)+P(K, I ) *BE(K, J) *WTG 
500 CONTINUE 

CALL MXINV(H, 18, IERH, NBETA) 

WRITE (6, 550)IERH 

550 FORMAT (/, ° Inversion Indicator Returned 
- By SSP(IBM) = °,I3) 

DO 625 1=1, NBETA 
DO 625 J=l,NDOF 
HG(I,J)=O.ODOO 
DO 600 K=1 , 1® ETA 

600 HG(I,J)=HG(I,J)+H(I,K)*G(K,J) 

625 HINVG(I,J)=HG(I,J) 

DO 675 1=1, NDOF 
DO 675 J=1,ND0F 
STMTX( I , J ) =0 . 0D00 
DO 650 K=l, NBETA 

650 STMTX(I,J)=STMTX(I,J)-K3(K,I)*HG(K,J) 

675 REALK(I,J)=STMTX(I,J) 

IF (10. EQ. 0) GO TO 800 
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WRITE (6,700 )NN0DES , NBETA 

700 FORMAT ( /, ° [ K ] Nfetrix Returned By § HISQUE t°/ 
SPECIFIED: Nutber of Nodes = °,I1,/ 

-° Number of Betas = 0 ,12) 

DO 725 KROW=l , NDOF 

725 WRITE(6, 775) (REALK(KRCW,KCOL) ,KC0Lf1 , NDOF) 

775 FORMAT ( / , IX, 8 ( El 3 . 7 , 3X) , / , IX, 8( El 3 . 7 , 3X) ) 

800 RETURN 


1000 ASTRNO “NODES 0 

BSTRNO=°OUT OF RANGE. 0 

WRITE ( 6 , 5000 ) ASTRNG, NN0DES , BSTRNG 

STOP 

2000 ASTRNG= “BETAS 0 

BSTRNG=“INSUFFICIENT. “ 

GO TO 4000 

3000 ASTRNG= “BETAS 0 

BSTRNG=“OUT OF RANGE. 0 
4000 WRITE (6, 5000 ) ASTRNG, NBETA, BSTRNG 
5000 FORMAT (IX, 

_«******************************************************* 0 t f 9 
-° * PROGRAM ABORT ACTIVATED EftT § HISQUE t : *“,/, 

-° * SPECIFIED NUMBER OF °,A5,° [°,I2,°] IS °,A13, ° *° 

J ' « *******************************************************° ) 

STOP 

END 
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SUBROUTINE QUEST (XS, YS, 15, 16, 17, 18,NBETA,HINVG,DISP, STRESS) 


DIMEM3ION DISP(l) , STRESS (1) ,HINVG(NBEJTA, 1) 
DIMENSION P(3, 18) ,PHINVG(3, 16) 

NN0DES=4+I5+I6+I7+I8 
NDOF =2*NNODES 

DO 100 IROW=l , 3 
DO 100 ICOLfI, 18 
100 P(IROW,ICOL)=O.ODOO 

P(1,1)=1.0D00 
P ( 1 , 4 ) =YS 
P (1 , 6)=XS 
P(1,8)=YS**2 
P(l, 10)=X5**2 
P(1,11)=2.D0*XS*YS 
P(1,13)=YS**3 
P(1,15)=XS**3 
P(1,16)=XS*(YS**2) 

P(1,17)=()S**2)*YS 

P(2, 2)=1. ODOO 

P ( 2 , 5 ) =XS 

P(2, 7)^ys 

P(2, 9)=>S**2 

P(2, 10)=Y5**2 

P(2,12)=2.D0*>5*YS 

P(2,14)=XS**3 

P(2,15)=3.D0*XS*(^S**2) 

P(2,17) = (YS**3)/3.D0 
P(2,18)=()S**2)*YS 

P(3, 3)=1 . ODOO 

P(3, 6)=-YS 

P(3,7)=-XS 

P ( 3 , 10 )=-2 . D0*XS *YS 

P(3,11)=-YS**2 

P(3,12)=-XS**2 

P(3,15)=-3.DO*(XS**2)*YS 

P(3,16)=(YS**3)/(-3.D0) 

P(3,17)=-XS*(YS**2) 

P(3,18)=()S**3)/(-3.D0) 

DO 200 1=1,3 
DO 200 J=1,ND0F 
PHINVG( I, J)=0. 0 
DO 200 K=1 , NBETA 
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PHINVG( I , J )■ =P ( I , K) *HINVG (K, J )+PHINVG( I , J ) 
200 CONTINUE 

DO 300 1=1,3 
STRESS ( I )= 0.0 

DO 300 K=1 , NDOF 

STRESS ( I ) =PHINVG ( I , K ) *DISP ( K ) +STRESS ( I ) 
300 CONTINUE 
RETURN 
END 
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B.2 


ELEMENT HEL710 


The following variables are required by subroutine HEL710: 

E =Young° s modulus 

NU =3?oisson° s ratio 

T = Element thickness 

X =Vector of nodal coordinates, x 

Y Sector of nodal coordinates, y 

lx =0(1) if node x is absent( present) 

10 Output control variable 

The following variables are returned by subroutine HEL710: 
HINVG =Product of [H]inv. and [G] matrices 
REALK =Element stiffness matrix 

Sibroutine H EI ST R carputes the stress vector STRESS from the 
following variables: 

XS,YS =Coordinates of point at Which stress is desired 
HINVG =Nbtrix returned by subroutine HEL710 

DISP ^Vector of nodal displacements 
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************** 


C KUROSH KAFIE MARCH 1981, M.I.T. * 

Q* ***** ************ ********* ************************************* 


SUBROUTINE HEL710(E,NU,T,X,Y,I8,I9,I10,IO,HINVG,REALK) 


IMPLICIT REAL*8 ( A-H » O-Z ) 

DOUBLE PRECISION NU,J11,J12,J21,J22 

REAL HINVG, REALK 

DIMENSION X(l) ,Y(1) , ZETA(5) , W(5 ) 

DIMENSION DIM(3 ) , F(10) ,DZ(10) ,DE(10) 

DIMENSION S(3,3),P(3,25),H(25,25) 

DIMEKEION BE(3, 20) ,G(25, 20) 

DIMENSION HG(25,20) ,SIMrX(20 / 20) 

DIMENSION HINVG(25, 1) ,REALK( (2* (7+I8+I9+I10) ),1) 

DDX(M)=+J22*DZ(M)-J12*DE(M) 

DDY(M)=-J21*DZ (M)+J11*DE(M) 

DATA ZETA/ -.9061 7984593 8664D0, 

2 -• 538469310105683D0, .00000000000000000, 

3 . 538469310105683D0, . 9061 7984593 8664D0/ 

DATA W/ . 2369268850561 89D0, 

2 . 4786286 70499366D0, . 56 888888888888 9D0, 

3 . 4 7862 86 704993 66D0, . 2369268850561 89D0/ 


D8 =18 
D9 =19 
D10=I10 

NNODES=7+I 8+1 9+110 
NDOF =2*NN0DES 

DO 150 IR=1, 25 

DO 100 ICH=1 , 2 5 
100 H(IR, ICH)=0. ODOO 
DO 125 ICG=1,20 
125 G (IR, ICG )=0. ODOO 
150 CONTINUE 

DO 250 IRC3W=1, 3 
DO 200 ICOL=l , 2 5 
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200 P ( I ROW , ICOL) =€) . 0D00 
DO 225 JCOLr=l ,20 
225 BE(IRCW,JCOL)=O.ODOO 
250 CONTINUE 

S(1,1)=1.D0/E 

S (1, 2)=-NU/E 

S(1,3)=O.ODO/E 

S (2,1 )=-NU/E 

S(2, 2)=1.D0/E 

S (2, 3)=0.0D0/E 

S(3,1)=O.ODO/E 

S(3,2)=0.0D0/E 

S ( 3 , 3 ) =2 . DO* (1 . D04NU) /E 

DO 500 11=1,5 
ZB=ZETA(II) 

DO 500 JJ=1, 5 
ET=ZETA( JJ) 

OPZ=l .DO-t-ZE 

OMZ=1.DO-ZE 

OPE=l . D04CT 

0ME>=1 . D0-E7T 

HPE=.5D0+ET 

HME=. 5D0-ET 

ESQ=d*ET 

ECU=ET*ESQ 

EQU=ET*ECU 

OMZSQ=1.DO-ZE*ZE 

CMESQ=1. DO-ESQ 

F( 8)= D8*OMZSQ*OME/2 . DO 
F( 9)= D9*CMESQ*0PZ/ 2 . DO 
F(10)=D10*a-1ZSQ*OPE/2 .DO 

F ( 1 ) =+HPE*CT*HME*CME*CMZ/3 . DO-F ( 8)/2.D0 
F(2)=OPZ*OME/4.DO-(F( 8)+F( 9))/2.D0 
F(3)=OPZ*OPE/4.DO-(F( 9)+F(10) )/2.D0 
F (4 ) =-OPE*HPE*ET*HME*CMZ/ 3 . DO-F (10) /2.D0 
F ( 5 ) =+OPE*HPE*ET *CME*CMZ/ . 7 5D0 
F ( 6 ) =+ 2 . DO*OPE*HPE*HME*OME*OMZ 
F( 7)=-0PE*ET*HME*CME*CMZ/ . 75DO 

XS=O.ODOO 

Y5=0.0D00 


DO 300 K=1 , NNODES 
XS=XS+X(K) *F(K) 
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YS=YS+Y(K)*F(K) 
300 CONTINUE 

XSQ=XS*XS 

XCU=)S*>5Q 

XQU=XS*XCU 

YSQ=YS*YS 

YCO=YS*YSQ 

YQU=YS *YCU 


P(1,1)=1.D00 

P(l, 4)=YS 

P(1,6)=XS 

P(1,8)=JSQ 

P(1,9)=2.D0*XS*YS 

P(1,10)=YSQ 

P(1,13)=XCU 

P(1,14)=(XSQ)*(YS) 

P(l» 15)=(XS )*(YSQ) 

P(1,16)=YCU 

P (1 , 19)=XQU 

P(1,20)=(XCU)*YS 

P(l» 21)=(XSQ) *(YSQ) 

P(1,22)=)S* (YCU) 

P(1,23)^YQU 


P(2, 2)=1.D00 

P(2,5)=XS 

P(2, 7)=YS 

P(2, 8)^YSQ 

P(2, 11)=2 .D0*XS*YS 

P(2, 12)=XSQ 

P(2,13)=3.DO*(XS)*(YSQ) 
P(2,14)=(YCU)/3.D0 
P ( 2 , 1 7 ) =XCU 
P(2,18)=(XSQ)*YS 
P(2, 19)=6. DO* (XSQ)* (YSQ) 
P(2, 20)=XS* (YCU) 
P(2,21)=(YQU) /6.0D00 
P(2, 24)=4. DO* (XCU) *YS 
P ( 2 , 2 5 ) =XQU 


P ( 3 , 3 ) =1 . DOO 

P(3,6)=-YS 

P(3 , 7)=-XS 

P(3, 8)=-2.D0*>S*YS 

P ( 3 , 9 ) =-YSQ 

P(3,11)=->5Q 

P(3,13)=-3.DO* (XSQ) * (YS ) 
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P(3,14)=-XS*(Y5Q) 

P(3,15)=-(YCU)/3.DO 
P(3,18)=-(XCU)/3.D0 
P(3,19)=-4.D0* (XCU ) *YS 
P(3,20)=-3.D0*(XSQ)*(Y5Q) /2.D0 
P ( 3 , 2 1 ) =-2 . DO*XS * (YCU ) /3 . DO 
P(3/22)=(YQU) /-4.D0 
P(3, 24)=-XQU 

DZ( 8)= -D8*ZE*CME 
DE( 8)= -D8*OMZSQ/2 . DO 
DZ( 9)= +D9*CMESQ/2 . DO 
DE( 9)= -D9*OPZ*ET 
DZ (10)=-D10*ZE*OPE 
DE(10)=+D1 0*CMZSQ/2 . DO 

DZ ( 1 ) = ( EQU-ECU- ( (ESQ-ET) /4.D0) )/(-3.D0)-DZ(8) /2.D0 
DE ( 1 ) = ( 4 . DO*ECU-3 . DO*ESCH-. 2 5DO-ET/2 . DO) *CMZ/3 .DO 
-DE(8)/2.D0 

DZ(2)=+CME/4.D0-0. 5DO*(DZ(8)+DZ (9) ) 
DE(2)=-OPZ/4.D0-0.5D0*(DE(8)+DE( 9)) 
DZ(3)=MOPE/4.D0-0.5D0*(DZ(9)+DZ(10) ) 
DE(3)=-KDPZ/4.D0-0.5D0*(DS(9)+DE(10) ) 

DZ (4)=(EQU+ECU-(ESQ+ET) /4.D0) / (-3 .DO )-DZ (10) /2 .DO 
DE(4) = (4. D0*ECLH-3.D0*ESQ-. 25D0-ET/2. DO) *OMZ/3.DO 
-DE(10) /2.D0 

DZ ( 5 ) = ( EQUf . 5DO*EO>-ESQ-ET/2 . DO ) / . 75D0 

DE ( 5 ) =- (4.DO*EO>1 . 5DO*ESQ-2 . DO*ET- . 5D0) *CMZ/ . 75DO 

DZ ( 6 ) =-2 . DO* ( EQU-1 . 2 5D0* ESCH- . 2 5D0 ) 

DE ( 6 )=+2 . DO* ( 4 . D0*ECU-2 . 5D0*ET ) *CMZ 

DZ(7)=+(BQU-.5D0*EQJ-ESQ+-ET/2.D0) /. 75D0 

DE ( 7 ) =- ( 4 . D0*ECT>1 . 5D0*ESQ-2 . DO*ET+ . 5D0) *CMZ/ . 7 5D0 

J11=0.0D00 
J12O.0D00 
J21=0. ODOO 
J22=0.0D00 

DO 350 K=1 , NNODES 
J11=DZ(K)*X(K)-KJ11 
J12=OZ(K) *Y(K)+J12 
J21=OE(K)*X(K)-hJ21 
J22=DE(K)*Y(K)-KJ22 
350 CONTINUE 

DETJ=vJ 1 1* J 22-J 1 2* J 2 1 

WTH=W( II ) *W (J J) *DETJ*T 
DO 400 J=l, 25 
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DO 370 K=l, 3 
DLM(K)=O.ODOO 
DO 370 L=l, 3 

370 DIM(K)=DLM(K)+S (K,L)*P(L, J) 

DO 390 I=J,25 
HTEMP=0 . 0D00 
DO 380 L=l, 3 

380 HTEMP=HTEWP+P(L,I)*DLM(L) 

390 H(I,J)=H(I,J)-WrEMP*WrH 
400 CONTINUE 

DO 435 1=1, (NDOF-1) , 2 

ITEMP=(I+l)/2 

BE(1,I)=DDX(ITEMP) 

BE(3,l)=DDY( ITEMP) 

435 CONTINUE 

DO 445 1=2 , NDOF , 2 
ITEKP=l/2 

BE(2,I)=DDY(ITEMP) 

BE(3,I)=DDX(ITEMP) 

445 CONTINUE 

WTG=W(II)*W(JJ)*T 
DO 450 1=1,25 
DO 450 J=l,NDOF 
EX) 450 K=l, 3 

450 G(I,J)=G(I,J)+P(K,I)*BE(K,J)*WTG 
500 CONTINUE 

CALL MXINV(H,25, IERH,0) 

WRTTE(6,550)IERH 

550 FORMAT (IX, ° INVERSION INDICATOR RETURNED 
- BY SSP(IBM) = 0 ,13) 

DO 625 1=1,25 
DO 625 J=l,NDOF 
HG(I,J)=O.ODOO 
DO 600 K=1 , 2 5 

600 HG( I, J)=HG(I, J)+H (I,K) ^ (K, J) 

625 HINVG ( I , J ) =HG ( I , J ) 

DO 675 1=1, NDOF 
DO 675 J=l,NDOF 
STMTX(I, J)=0.0D00 
DO 650 K=1 , 2 5 

650 STMTX(I, J)=STMTX(I,J)-K5(K, I)*HG(K, J) 
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675 REALK(I,J)=STMTX(I,J) 

IF(IO.BQ.O) GO TO 999 
WRTTE(6,900) 

DO 700 KR=1 , NDOF 

700 WRTTE(6,800) (REALK(KR, KC) ,KC=l,NDOF) 

800 FORMAT ( / , ° \10E13. 5,/, ° MQE13.5) 

900 FORMAT ( / / , ° [ K ] matrix returned by § HEL710 t *//) 


999 RETURN 
END 
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SUBROUTINE HELSTR(XS , YS , 18, 19 , 110, HINVG, DISP, STRESS) 


DIMENSION P(3,25) ,HINVG(25, 1) ,PHINVG(3,20) 
DIMENSION DISP(2* (7+I8+I9+I10) ) ,STRESS(3) 

DO 100 IROW=l, 3 
DO 100 ICOLfI, 25 
100 P(IROW,ICOL)=O.DO 


P(l, 1)=1. D00 
P(1,4)=YS 
P ( 1 , 6 ) =XS 
P(1,8)=X3**2 
P(1,9)=2.D0*XS*YS 
P(1,10)=YS**2 
P(l, 13)=XS**3 
P(1,14)=(XS**2)*(YS) 

P(1 , 15) = (XS ) *(YS**2 ) 

P(1,16)=^S**3 

P(l, 19)=XS**4 

P(1,20)=(XS**3)*YS 

P(1 , 21) = (XS**2 ) *(YS**2 ) 

P(1,22)=)S* (YS**3) 

P(l, 23)=YS**4 


P(2,2)=1.D00 

P(2, 5)=XS 

P ( 2 , 7 ) =YS 

P(2, 8)=YS**2 

P(2,ll)=2.DO*XS*YS 

P(2, 12)=XS**2 

P ( 2 , 1 3 ) =3 . DO* ( >5 ) * ( YS* *2 ) 

P(2,14) = (YS**3)/3.D0 

P(2,17)=XS**3 

P(2, 18) = (XS**2 ) *YS 

P(2,19)=6.D0*(XS**2)*(YS**2) 

P(2, 20)=XS* (YS**3 ) 

P(2,21)=(YS**4) /6.0D00 

P(2,24)=4.D0*(XS**3)*YS 

P(2,25)=>S**4 


P(3, 3)=1.D00 

P ( 3 , 6 ) =-YS 

P(3, 7)=-XS 

P(3,8)=-2.D0*>5*YS 

P(3,9)=-YS**2 

P(3,11)=-XS**2 

P(3, 13)=-3.D0* (XS**2 ) * (YS ) 

P(3,14)=-XS*(YS**2) 
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P(3,15)=-(YS**3)/3.DO 
P(3,18)=-(XS**3)/3.D0 
P (3 , 19)=-4. DO* (XS**3 ) *YS 
P(3,20)=-3.D0*(X3**2)*(YS**2) /2.DO 
P ( 3 , 2 1 ) =-2 . DO *XS * ( YS * *3 ) /3 . DO 
P(3,22)=(YS**4)/-4.D0 
P (3 , 24) =-l . DO* (XS **4) 

NNODES=7+I8+ 19+11 0 
NDOF =2*NNODES 

DO 200 1=1,3 
DO 200 J=l, NDOF 
PHINVG(I,J)=0.0 
DO 200 K=1 , 2 5 

PHINVG ( I , J ) =P ( I , K ) *H INVG ( K, J ) +PHINVG ( I , J ) 
200 CONTINUE 

DO 300 1=1,3 
STRESS (I )=0.0 
DO 300 K=1 , NDOF 

STRESS ( I ) =PHINVG( I , K ) *DISP ( K ) +STRESS ( I ) 
300 CONTINUE 
RETURN 
END 
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3.3 ELEMENT PET4/12 


The following variables are required by subroutine PET4/12: 

E =Young°s modulus 

NU =€>oisson°s ratio 

T = Element thickness 

X ^Vector of nodal coordinates, x 

Y =Vector of nodal coordinates, y 

10 Output control variable 

The stiffness matrix, REALK, is returned by subroutine PETT4. 
Subroutine PETSAP ccnputes the stress vector STRESS frcm the 
following variables: 

XS,YS Coordinates of point at which stress is desired 
RZ = Radius of the circular hole 

DISP =Vector of nodal displacements 
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c********************************************************^*^^ 

C KURCSH KAFIE DECEMBER 1980, M.I.T.* 

C*************************************************************** 


SUBROUTINE PET4(E,NU,T,X, Y,IO,REALK) 

C*************************************************************** 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 


DESCRIPTION: * 

SUBROUTINE PET4 DETERMINES THE STIFFNESS MATRIX OF A * 

FOUR NODE MEMBRANE ELEMENT WITH A BUILT IN CIRCULAR * 

TRACTION-FREE EDGE: * 

* 

E =YOUNG°S MODULUS. * 

NU=PO ISSON 0 S RATIO. 0 3 * 

T =T HICKNESS. . . * 

. . * 

. . * 

. . * 

4 0. * 

• . * 

. . ^ 

o. . . o * 

12 * 

* 


C*********^***************************************************** 


IMPLICIT REAL*8( A-H, 0-Z ) 


REAL HINVG, REAUC(8, 1 ) 

COMMON /PET1X2/HINVG(12,8) 

DOUBLE PRECIS ION J, LNGT , L, NU, NUX, NUY, NUT 

DIMENSION X(1),Y(1),GS(7),GW(7),DIM(3) 
DIMENSION S (3, 3) ,P(3, 12) ,H(12, 12) 
DIMENSION NUT (3, 2) ,L(2, 8) ,G(12, 8) 
DIMENSION HG(12, 8) ,STMIX(8,8) 


DATA GS/ 

2 -0. 741531185599394D0, 

3 0 . 00000000000000 0D0 , 

4 0. 741531185599394D0, 

CftTA GW/ 

2 0.2 797053914892 77DO, 

3 0.417959183673469D0, 

4 0. 2797053914892 77DO, 

NDOF = 8 
NBETA=12 


-0.9491079123 42759D0, 
-0. 405845151377397D0, 
0. 40584515137739 7D0, 
0. 94910 791 2342 759DO/ 

0 . 1 294849661 68870D0 , 
0.381 83 0050505119D0, 
0 . 3 81 830050 5051 1 9D0 , 
0. 1 294849661 6887QD0/ 


DO 75 IR=1 , NBETA 
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DO 25 ICH=1, NBETA 
25 H(IR, ICH)=O.DOO 
DO 50 103=1, NDOF 
50 G(IR,ICG)=0.D00 
75 CONTINUE 

R2HDSQRT(X(1 ) **2+Y(l ) **2 ) 

CALL ANGLE(X(2),Y(2),THETA2) 
CALL ANGLE(X(3),Y(3) ,THETA3) 
THETA0= (THETA3 -THETA2 ) /2 . ODOO 
ALPHA =(THETA3-KTHETA2)/2.0D00 
C1=(X(2)*Y(3)-X(3)*Y(2) ) 
C2=X(2 )-X(3 ) 

C3=Y(2)-Y(3) 

IFLAG=0 

S(1,1)=1.D0/E 

S(l, 2)=-NU/E 

S (1, 3)=0.0D0 

S (2, 1 )=-NU/E 

S(2,2)=1.D0/E 

S(2,3)=O.ODO 

S ( 3 , 1 ) =0 . ODO 

S(3,2)=0.0D0 

S (3, 3 )=2 .DO* (1 .DO+NU) /E 


DO 300 11=1,7 
ZE=GS (II ) 

DO 300 JJ=1, 7 
ET=GS (JJ) 

THETA=ZE*THETAO+ALPHA 

PHI=Cl/(C2*DSIN(THETA)-C3*DCOS (THETA) ) 

R=( (PHI-RZ) /2 . ODO) *ET+ (PHI+RZ) /2 . ODO 
J=THETAO*( ( (PHI-RZ) **2) *ET+PHI**2-RZ**2) /4.D0 

100 ST=DS IN (THETA) 

CT =DCOS (THETA) 

P(1,1)=2.D0*(1.D0-RZ/R) 

P(l, 2)=(1.D0-(RZ**2/R**2) )*CT 
P(1,3) = (1.D0-(RZ**2/R**2) )*ST 
P(1,4)=3.D0*(R-(RZ**2/R) ) 

P(l, 5)=2.D0*(R-(RZ**3/R**2) )*CT 
P(1,6)=2.D0*(R-(RZ**3/R**2) )*ST 
P(1 , 7) =4. DO* (R**2-( RZ**3/R) ) 

P ( 1 , 3 ) =3 . DO* ( R* *2- ( RZ* *4/R* *2 )) *CT 
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P(1,9)=3.D0*(R**2-(RZ**4/R**2) )*ST 
P( 1,10) =5 .DO* (R**3-(RZ**4/R) ) 

P ( 1 , 1 1 ) =4 . DO* ( R* *3 -( RZ* *5 / R* *2 )) *CT 
P (1 ,12)=4.D0* (R**3-(RZ**5/R**2) ) *ST 

P ( 2 , 1 ) =2 . DO 
P(2,2)=2.D0*CT 
P(2, 3)=2.D0*ST 
P (2, 4)=6.D0*R 
P(2,5)=6.D0*R*CT 
P(2,6)=6.D0*R*ST 
P(2, 7) =12. DO* (R**2 ) 

P (2, 8) =12. DO* (R**2)*CT 
P(2,9)=12.D0*(R**2)*ST 
P (2, 10) =20. DO* (R**3 ) 
P(2,11)=20.D0*(R**3) *CT 
P( 2, 12) =20. DO* (R**3) *ST 

P(3,1)=0.D00 

P(3,2)=(1.D0-(RZ**2/R**2))*ST 
P(3,3)=( (RZ**2/R**2)-1. DO) *CT 
P(3,4)=O.DOO 

P ( 3 , 5 ) =2 . DO* (R- ( RZ* *3/R**2 ) ) *ST 
P ( 3 , 6 ) =2 . DO* (( RZ* *3/R* *2 ) -R) *CT 
P(3, 7)=0.D00 

P ( 3 , 8) =3 . DO* (R**2-( RZ* *4/R* *2 ) ) *ST 
P(3,9)=3.DO*( ( RZ* *4/R* *2 ) -R* *2 ) *CT 
P(3,10)=0.D00 

P ( 3 , 1 1 ) =4. DO* ( R* *3 - ( RZ* *5/R**2 ) ) *ST 
P(3, 12)=4.D0*( (RZ**5/R**2)-R**3)*CT 

IF(IFLAG.EQ.I) GOTO 400 

WTH=J*OV(ll)*GW(JJ) 

DO 200 KC=1 , NBETA 
DO 125 NN=1, 3 
DLM(NN)=O.DOO 
DO 125 KR=1,3 

125 DLM(NN)=DLM(NN)+S (NN,KR)*P(KR,KC) 

DO 175 KR=KC,NBETA 
HTEMP=O.ODOO 
DO 150 NN=1,3 

150 HTEMP=OTE>P+P(NN,KR)*DIM(NN) 

175 H(KR,KC)=«(KR,KC)4HrEMP*WrH 
200 CONTINUE 

300 CONTINUE 
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CALL MXINV(H,NBETA, IERH, 0) 


IFLAG=1 

DO 800 IX=1,3 
XA=X( IX) 

XB=X(IX+1) 

YA=Y ( IX) 

YB=¥ ( IX+1 ) 

LNGT=DSQRT ( (XA-XB ) **2+ (YA-YB ) **2 ) 
CALL ANCLE ( (YB-YA) , (XA-XB) ,CMEGA) 
NUX=DCOS (OMEGA) 

NUY=D6IN(CMEGA) 


DO 700 11=1,7 

XS=GS (II)* (XB-XA ) /2 . 0D0+(XB+XA ) /2 . 0D0 
YS=GS ( II ) * ( YB-YA) /2 . 0D0+ ( YB+YA) /2 . ODO 

350 R=DSQRT (XS**2+YS**2 ) 

CALL ANGLE(XS,YS, THETA) 


GO TO 100 


400 NUT(1,1)=NUX*(CT**2)+NUY*CT*ST 
NUT ( 2 , 1 ) =NUX* (ST* *2 ) -NUY*CT*ST 
NUT(3, 1 )=3sTUY* (CT**2-ST**2 )-2.D0*NUX*CT*ST 

NUT ( 1 , 2 ) =NUY* (ST* *2 ) 4NUX*CT*ST 
NUT(2, 2)=£IUY* (CT**2 )-NUX*CT*ST 
NUr(3,2)=NUX*(CT**2-6T**2)+2.D0*NUY*CT*ST 

DO 500 I ROW=l , 2 
DO 500 ICOL=l , 8 
L ( I RCW , ICOL) =0 . 0D00 
500 CONTINUE 


ENT1= (1 . DO-GS ( 1 1 ) ) /2 . ODO 
ENT2= ( 1 . D0-K3S ( 1 1 ) ) /2 . ODO 

L(l, (2*IX-1) )=ENT1 
L(l, (2*IX+1) )=ENT2 

L(2, (2*IX) ) =ENT1 


129 


L(2, (2*IX+2) )=ENT2 

WTG=GW ( 1 1 ) *LNGT/2 . DO 

DO 600 KC=1,ND0F 
DO 525 NN=1,3 
DOM ( NN ) =0 • D00 
DO 525 KR=1, 2 

525 DOM(NN)=C>OM(NN)+NUr(NN,KR)*L(KR,KC) 

DO 575 KR=1 , NBETA 
GTEMP=0. ODOO 
DO 550 NN=1, 3 

550 GTEMP=GTEMP+P(NN,KR)*DLM(NN) 

575 G(KR, KC)=G(KR,KC )-KjTEMP*WTG 
600 CONTINUE 

700 CONTINUE 

800 CONTINUE 

DO 850 M=1 , NBETA 
DO 850 N=l,NDOF 
HG(M,N)=0. ODOO 
DO 825 K=l, NBETA 

825 HG(M,N)=HG(M,N)+H(M, K) *G(K,N) 

850 HINVG(M,N)=HG(M,N) 

DO 900 M=l,NDOF 
DO 900 N=l, NDOF 
STMTX(M,N)=0.0D00 
DO 875 K=l, NBETA 

875 STMTX(M / N)=STMTX(M > N)+G(K,M)*HG(K,N) 

900 REALK(M / N)=STMTX(M,N) 

IF(IO.EQ.O) GO TO 1000 
WRITER, 975) 

DO 925 KR=1 , 8 

925 WRITE(6, 950) (REALK(KR,KC) ,KC=1, 8) 

950 FORMAT (IX, 8( El 4. 5, 2X)//) 

975 FORMAT ( //, ° [ K ] Matrix Returned By Subroutine § PET4/12 t 

V/) 

1000 RETURN 
END 
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SUBROUTINE PETSAP(XS,YS, RZ,DISP, STRESS) 
COMMON /PET1X2/HINVG (12,8) 

DIMENSION DISP(l) ,STRESS(1) 

DIMENSION P ( 3 , 1 2 ) , PHINVG (3,8) 

DOUBLE PRECISION XD,YD,THETAD 

R=SQRT (XS**2+YS**2) 

XD=XS 

YD=YS 

CALL ANGLE(XD,YD,THETAD) 

THETA=THETAD 

ST=SIN(THETA) 

CT=COS (THETA) 

P ( 1 , 1 ) =2 . 0* (1 . 0-RZ/ R) 
P(1,2)=(1.0-(RZ**2/R**2) )*CT 
P(1,3)=(1.0-(RZ**2/R**2) )*ST 
P(1,4)=3.0*(R-(RZ**2/R) ) 

P(1 , 5)=2. 0* (R-(RZ**3/R**2) )*CT 
P(1,6)=2.0*(R-(RZ**3/R**2))*ST 
P(l f 7)=4.0*(R**2-(RZ**3/R) ) 
P(1,8)=3.0*(R**2-(RZ**4/R**2) )*CT 
P(1,9)=3.0*(R**2-(RZ**4/R**2))*ST 
P(1,10)=5.0*(R**3-(RZ**4/R) ) 

P ( 1 , 1 1 ) =4. 0* ( R**3 -( RZ* *5 /R**2 )) *CT 
P ( 1 , 12 ) =4. 0* ( R* *3- ( RZ* *5/R* *2 ) ) *ST 


P (2,1) =2.0 
P ( 2 , 2 ) =2 . 0*CT 
P(2, 3)=2.0*ST 
P (2, 4)=6.0*R 
P(2, 5)=6.0*R*CT 
P (2, 6)=6.0*R*ST 
P(2, 7)=12.0*(R**2) 

P(2, 8)=12.0*(R**2)*CT 
P(2, 9)=12. 0* (R**2 ) *ST 
P ( 2 , 10 ) =2 0 . 0* ( R* *3 ) 
P(2,11)=P(2,10)*CT 
P(2,12)=P(2,10)*ST 

P(3, 1)0.00 

P(3,2)=(1.0-(RZ**2/R**2) )*ST 
P(3/3) = ( (RZ**2/R**2)-1.0)*CT 
P(3, 4)=0.00 

P(3 , 5)=2. 0* (R— ( RZ**3/R**2 ) )*ST 
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P(3,6)=2.0*( (RZ**3/R**2 )-R) *CT 
P(3,7)=0.00 

P ( 3 , 8) =3 . 0* ( R* *2 -( RZ* *4/R**2 ) ) *ST 
P(3,9)=3.0*( (RZ**4/R**2 )-R**2) *CT 
P(3,10)=0.0 

P(3, 11)=4.0*(R**3-(RZ**5/R**2) )*ST 
P(3,12)=4.0*( (RZ**5/R**2)-R**3)*CT 

DO 100 1=1,3 
DO 100 J=l, 8 
PHINVG( I, J)=0. 0 
DO 100 K=l,12 

PHINVG( I , J ) =P ( I , K ) *HINVG(K, J )+PHINVG( I , J ) 
100 CONTINUE 

DO 200 1=1,3 
STRESS (I )=0,0 
DO 200 K=l, 8 

STRESS ( I ) =PHINVG( I , K) *DISP (K)-tSTRESS ( I ) 
200 CONTINUE 

RETURN 

END 
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B.4 ELEMENT PET4X5 


The following variables are required by subroutine Pr7T4X5: 

E =Young° s modulus 

NU Poisson® s ratio 

T = Element thickness 

X =Vector of nodal coordinates, x 

Y =Vector of nodal coordinates, y 

10 =Output control variable 

The stiffness matrix, REftLK, is returned by subroutine PET 4. 
Subroutine PETSAP ccrrputes the stress vector STRESS from the 
following variables: 

XS,YS =Coordinates of point at which stress is desired 
RZ =Radius of the circular hole 

DISP Sector of nodal displacervents 


133 


C***************************** ********************************** 

C KURDSH KAFIE DECEMBER 1980, M.I.T.* 

C*************************************************************** 


SUBROUTINE PE7T4(E,NU,T,X, Y, IO,REALK) 


IMPLICIT REAL*8( A-H,0-Z) 

REAL HINVG, REALK (8,1) 

COMMON /PE7T1X2/HINVG(5, 8) 

DOUBLE PRECISION J,LNGT, L,NU,NUX,NUY,NUT 

DIMENSION X(1),Y(1),GS(7),GW(7),DIM(3) 
DIMENSION S(3,3),P(3,5),H(5,5) 

DIMENSION NUT(3 , 2) ,L(2, 8) ,G(5, 8) 

DIMENSION HG(5,8),STMTX(8,8) 

DATA GS/ -0. 949107912342759D0, 

2 -0. 741531185599394D0, -0. 405 8451 51 3 7 73 9 7D0, 

3 0.00000000000000000, 

4 0. 741531185599394D0, 

DATA GW/ 

2 0. 2797053914892 77DO, 

3 0.41 79591 8367346SD0, 

4 0. 279705391 4892 7 7D0, 

NDOF =8 
NBETA=5 

DO 75 IR=1, NBETA 

DO 25 ICH=1,NBETA 
25 H(IR,ICH)=O.DOO 
CO 50 I 03=1, NDOF 
50 G(IR,ICG)=0.D00 
75 CONTINUE 

RZ=OSQRT(X(l ) **2+Y(l ) **2 ) 

CALL ANGLE (X (2) ,Y(2) ,THEJTA2) 

CALL ANGLE(X(3),Y(3),THETA3) 

THEJTA0= (THET A3 -THETA2 ) /2 . 0D00 
ALPHA =(THETA3-ETHETA2)/2.0D00 
C1=(X(2)*Y(3)-X(3)*Y(2) ) 

C2=X(2 )-X(3 ) 

C3=Y(2 )-Y(3 ) 

IFLAG=0 


0. 405845151377397D0, 
0 . 94910 7912342 75 9D0/ 

0. 129484966168870D0, 
0 . 3 81 83005050511 9D0 , 
0. 3 81 83005050511SD0, 
0 . 129484966168870D0/ 
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S(1,1)=1.D0/E 

S(l, 2)=-NU/E 

S(1,3)=O.ODO 

S (2, 1)=-NU/E 

S (2, 2)=1.D0/E 

S (2 , 3)=0. ODO 

S ( 3 , 1 ) =0 • ODO 

S(3,2)=O.ODO 

S (3, 3) =2. DO* (1 .DO-HSTU) /E 


DO 300 11=1, 7 
ZE=GS(II) 

DO 300 JJ=1, 7 
ET=GS (JJ) 

THETA=ZE*THETAD+ALPHA 

PHI=C1/(C2*D6IN(THETA)-C3*DC0S (THETA) ) 

R= ( (PHI-RZ) /2 . ODO ) *ET+ (PHI+RZ) /2 . ODO 
J=THETAO*( ( (PHI-RZ)**2) *ET+PHI**2-RZ**2) / 4. DO 

100 ST=DSIN(THETA) 

CT =DCOS (THETA) 

S2T=OSIN(2.DO*THETA) 

C2T=DC0S (2.D0*THETA) 

P(1,1)=(1.D0-(RZ**2/R**2) ) 

P(1,2)=(R-(RZ**4/R**3) )*CT 
P(1,3)=(R-(RZ**4/R**3) )*ST 

P ( 1 , 4)= ( 1 • DO+3 . DO* ( RZ* *4/R* *4 )-4 . DO* ( RZ* *2/R* *2) ) *C2T 
P ( 1 , 5 ) = (1 . DO+3. DO* (RZ* *4/R**4 )-4. DO* (RZ**2/R**2 ) ) *S2T 

P(2, 1)=(1.D0+(RZ**2/R**2) ) 

P(2,2)=(3.D0*Rf(RZ**4/R**3) )*CT 
P(2, 3) = (3.D0*RH-(RZ**4/R**3) )*ST 
P(2,4)=(-1.D0-3.D0*(RZ**4/R**4) )*C2T 
P(2, 5)=(-l.D0-3.D0* (RZ**4/R**4) ) *S2T 

P(3, 1)=O.ODOO 
P(3,2)=(R-(RZ**4/R**3))*ST 
P(3, 3) = ( (RZ**4/R**3)-R) *CT 

p ( 3 , 4)= (-1 . DO+3 . DO* ( RZ* *4/R**4 )-2 . DO* (RZ* *2/R**2) ) *S2T 
P ( 3 , 5 ) = (+1 • DO-3 . DO* ( RZ* *4/R* *4 ) +2 . DO* ( RZ* *2 /R**2 ) ) *C2T 

IF ( IFLAG. EQ . 1 ) GO TO 400 

WTH=J*GW(II)*GW(JJ) 

DO 200 KC=1,I©ETA 
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DO 125 NN=1,3 
DLM(NN)=O.DOO 
DO 125 KR=1,3 

125 DLM(NN)=DLM(NN)+S (NN,KR) *P (KR,KC ) 
DO 175 KR=KC,NBETA 
HTEMP=0 . 0D00 
DO 150 NN=1 , 3 

150 ffTEMP=HTn<P+P ( NN, KR ) *DLM (NN) 

175 H(KR,KC)=H(KR,KC)4HTEMP*WTH 
200 CONTINUE 

300 CONTINUE 

CALL MXINV(H,NBETA,IERH,0) 

IFLAG=1 

DO 800 IX=1 , 3 
XA=X(IX) 

XB=X( IX+1 ) 

YA=Y( IX) 

YB=¥(IX+1) 

LNGT=DSQRT ( (XA-XB ) **2+ (YA-YB ) **2 ) 
CALL ANGLE ( (YB-YA), (XA-XB) ,CMEGA) 
NUX=C)COS( OMEGA) 

NUY=DSIN(CMEGA) 


DO 700 11=1, 7 

XS=GS (II)* (XB-XA ) /2 . 0D0+(XB+XA ) /2 . 0D0 
YS=GS ( II ) * ( YB-YA) /2 . 0D0+ ( YB+YA) /2 . 0D0 

350 R=DSQRT(XS**2+YS**2) 

CALL ANGLE (XS,YS, THETA) 


GO TO 100 


400 NUT(1, 1)=NUX*(CT**2)4NUY*CT*ST 
NUT(2,1)^UX*(ST**2)-NUY*CT*ST 
NUT(3, 1)=NUY* (CT**2-ST**2)-2.D0*NUX*CT*ST 

NUT (1,2) ^JUY* (ST* *2 ) +NUX*CT*ST 
NUT (2 , 2 ) =NUY* ( CT**2 ) -NUX*CT*ST 
NUT ( 3 , 2 ) =NUX* ( CT * *2 -ST ** 2 ) +2 . D0*NUY*CT *ST 
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DO 500 IR0W=1, 2 
DO 500 ICOL=l, 8 
L ( IROW, ICOL)=0 . 0D00 
500 CONTINUE 


ENT1=(1 .DO-GS (II) )/2.0D0 
ENT2= ( 1 . D0-K3S ( 1 1 ) ) /2 . ODO 

L(l, (2*IX-1) )=ENT1 
L(l, (2*IX+1))=ENT2 

L(2, (2*IX) ) =ENT1 
L(2, (2*IX+2) )=ENT2 

WTOGW ( 1 1 ) *LNGT/2 . DO 

DO 600 KC=l,NDOF 
DO 525 NN=1, 3 
dlm(nn)=o.doo 
DO 525 KR=1, 2 

525 DlM(NN)=CllM(NN)+NUr(NN,KR)*L(KR,KC) 

DO 575 KR=1 , NBETA 
GTEMP=0.0D00 
DO 550 NN=1,3 

550 GTEMP=<jTEMEH-P (NN,KR) *DLM(NN) 

575 G(KR,KC)=G(KR,l<C)-K7rEMP*WrG 
600 CONTINUE 


700 CONTINUE 

800 CONTINUE 

DO 850 M=l, NBETA 
DO 850 N=l, NDOF 
HG(M,N)=0.0D00 
DO 825 K=l, NBETA 

825 HG(M,N)=tiG(M,N)+H(M,K) *G(K, N) 

850 HTNVG(M,N)=HG(M,N) 

DO 900 M=1 , NDOF 
DO 900 N=l, NDOF 
STMTX(M,N)=O.ODOO 
DO 875 K=l, NBETA 

875 S'IMrX(M,N)=6TMTX(M,N)+G(K, M) *HG(K f N) 
900 REALK(M,N)=STMTX(M,N) 
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IF(IO.EQ.O) GO TO 1000 
WRTTE(6 f 975) 

DO 925 KR=1, 8 

925 WRTTE(6, 950) (REALK(KR,KC) ,KC=1, 8) 

950 F0RMAT(1X,8(E14.5,2X)//) 

975 FORMAT (//, ° [ K ] Matrix Returned By Subroutine § PET4X5 t 

7 /) 

1000 RETURN 
END 
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SUBROUTINE PETSAP(XS,YS, RZ,DISP, STRESS) 

CCMMDN /PET1X2/HINVG (5, 8) 

DIMENSION DISP{1) ,STRESS(1) 

DIMENSION P(3, 5) ,PHINVG(3, 8) 

DOUBLE PRECISION XD, YD,THETAD 

R=SCFT(>E**2+YS**2) 

XD=XS 

YD=YS 

CALL ANGLE(XD, YD, THETAD) 

THETA=THETAD 

ST=SIN(THETA) 

CT=COS (THETA) 

S2T=DSIN(2 . DO*THETA) 

C2T=DC0S (2.D0*THETA) 

P(1,1)=(1.0-(RZ**2/R**2)) 

P(1,2)=(R-(RZ**4/R**3) )*CT 
P (1 , 3 ) = ( R- ( RZ**4/R**3 ) ) *ST 

P ( 1 , 4) = ( 1 .0+3 . 0* ( RZ* *4/R* *4 )-4. 0* ( RZ* *2/R* *2) ) *C2T 
P (1 , 5)=(1 . 0+3. 0* (RZ**4/R**4)-4. 0* (RZ**2/R**2 ) ) *S2T 

P(2, 1)=(1.0+(RZ**2/R**2) ) 

P{2, 2)=(3. 0*R+(RZ**4/R**3 ) )*CT 
P(2,3)=(3.0*R+(RZ**4/R**3))*ST 
P(2, 4)=(-l. 0-3. 0* (RZ**4/R**4) )*C2T 
P(2 f 5)=(-1.0-3.0*(RZ**4/R**4))*S2T 

P(3, 1)=0.000 

P (3, 2)=(R-(RZ**4/R**3) ) *ST 
P(3, 3)=( (RZ**4/R**3)-R )*CT 

P ( 3 , 4 )= (-1 .0+3 . 0* ( RZ* *4/R* *4 )-2 . 0* ( RZ* *2/R* *2 ) ) *S2T 
P(3, 5)=(+l. 0-3. 0*(RZ**4/R**4)+2. 0*(RZ**2/R**2) )*C2T 


DO 100 1=1,3 
DO 100 J=l, 8 
PHINVG ( I , J ) =0 . 0 
DO 100 K=l, 5 

PHINVG( I, J ) =P ( I , K ) *HINVG(K, J )+PHINVG( I , J) 
100 CONTINUE 


DO 200 1=1,3 
STRESS (I )=0.0 
DO 200 K=l, 8 
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STRESS ( I ) =PHINVG ( I , K ) *DISP(K)+STRESS ( I ) 
200 CONTINUE 

RETURN 

END 
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